heun's method in matlab

main.m
function [t, y] = heun(dydt, tspan, y0, h)
% Implements Heun's method for solving ordinary differential equations
% dydt: function handle of the differential equation (in the form: dydt = @(t,y) ...)
% tspan: [t0, tf] where t0 is the initial value of t and tf is the final value of t
% y0: initial value of y given at t0
% h: step size
% t: vector of times at which the solution was computed
% y: vector of the solution at each corresponding time in t

t0 = tspan(1);
tf = tspan(2);
t = t0:h:tf;
y = zeros(1, length(t)); % preallocate solution vector
y(1) = y0;

for i = 1:length(t)-1
    k1 = dydt(t(i), y(i)); % evaluate slope at current point
    k2 = dydt(t(i)+h, y(i)+h*k1); % evaluate slope at point extrapolated from initial point using k1
    y(i+1) = y(i) + (h/2)*(k1 + k2); % update solution using trapezoidal rule
end
829 chars
21 lines

main.m
dydt = @(t,y) y - t^2 + 1; % define differential equation
tspan = [0, 2]; % specify time interval
y0 = 0.5; % specify initial condition
h = 0.25; % specify step size
[t, y] = heun(dydt, tspan, y0, h); % call heun function
plot(t,y) % plot solution
248 chars
7 lines