function [t,y] = rungekutta4( f, tspan, y0, n)
%RUNGEKUTTA4 Apply Runge-Kutta order 4 method to solve y' = f(y,t), y(a) = y0,
%  on the interval tspan(1) <= t <= tspan(2) with n steps.
%
%  Example usage:
%    f = @(t,y) t.*y.^2;
%    [t,y] = rungekutta4(f,[0,1],1,20);
%    plot(t,y);
%
    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
        m1 = f(t(j), y(j));  % slope at current point
        m2 = f(t(j)+h/2, y(j)+(h/2)*m1);  % slope at (approx) midway point
        m3 = f(t(j)+h/2, y(j)+(h/2)*m2);  % slope at (better approx) midway point
        m4 = f(t(j+1), y(j)+h*m3);  % slope at (approx) next point
        y(j+1) = y(j) + h*(m1+2*m2+2*m3+m4)/6;  % move by weighted average of slopes
    end
end