midpoint method in matlab

main.m
function [t, y] = midpoint(f, tspan, y0, h)
% Inputs:
%   f     - function handle for the right-hand side of the ODE
%   tspan - [t0, tf], where t0 is the initial time and tf is the final time
%   y0    - initial value for the solution
%   h     - step size for the numerical method
% Outputs:
%   t     - vector of time steps
%   y     - vector of corresponding solution values

t0 = tspan(1);
tf = tspan(2);
t = t0:h:tf;
y = zeros(size(t));
y(1) = y0;

for i = 1:(length(t)-1)
    k1 = f(t(i), y(i));
    k2 = f(t(i) + h/2, y(i) + h/2*k1);
    y(i+1) = y(i) + h*k2;
end

end
577 chars
24 lines