function [t,y] = midpoint( f, tspan, y0, n)
%MIDPOINT Apply the midpoint method to solve y' = f(y,t), y(a) = y0,
%  on the interval a <= t <= b with n steps.
%  This is a second order Runge-Kutta method, and is similar to the
%  midpoint rule for quadrature.
    a = tspan(1);
    b = tspan(2);
    h = (b-a)/n;
    t = a:h:b;
    y = zeros(1,n+1);
    y(1) = y0;
    for j = 1:n
        F1 = f(t(j),y(j));  % slope at current point
        F2 = f(t(j)+.5*h,y(j)+.5*h*F1);  % slope at approx. mid point
        y(j+1) = y(j) + h*F2;  % Move by midpoint slope.
    end
end