function y = rk2_1d( 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(0) = y0 and time t = t0.  return the entire vector
%  y = [ y(0), y(1), ..., y(n_steps) ]
%

y = zeros( n_steps+1, 1);  %  create an empty vector to put answers in

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
