function y = rk2( y0, t0, dt, n_steps, deriv_func ) % % take n_steps time steps, of size dt, in the differential equation % dy/dt = deriv_func( y, t ) % begining with y(1) = y0 and time t = t0. % all quantities yand f are of dimension m. m = length( y0 ); % detemine dimension from y0 y = zeros( n_steps+1, m ); % create an empty array for results y( 1, : ) = y0; % initial condition for n = 1:n_steps % the first half-step t = t0 + dt*(n-1); f = feval( deriv_func, y(n,:), t ); y_star = y(n,:) + 0.5*dt*f; % the full step t_star = t + 0.5*dt; f_star = feval( deriv_func, y_star, t_star ); y(n+1,:) = y(n,:) + dt*f_star; end