function [t, u] = ab4(f, tspan, u0, n)
% AB4     4th-order Adams-Bashforth formula for an IVP.
%
% Input:
%   f       defines f in u'(t)=f(t,u) (function)
%   tspan   endpoints of time interval (2-vector)
%   u0      initial value (m-vector)
%   n       number of time steps (integer)
% Output:
%   t       selected nodes (vector, length n+1)
%   u       solution values (array, size n+1 by m)
%
    % Discretize time.
    a = tspan(1);  
    b = tspan(2);
    h = (b-a)/n;
    t = tspan(1) + (0:n)'*h;

    % Constants in the AB4 method.
    k = 4;  
    sigma = [55; -59; 37; -9]/24;  

    % Find starting values by RK4.
    u = zeros(length(u0), n+1);
    [ts, us] = rk4(f,[a, a+(k-1)*h],u0,k-1);
    u(:,1:k) = us(1:k,:).';

    % Compute history of u' values, from oldest to newest.
    ff = zeros(length(u0),k);
    for i = 1:k-1
        ff(:,k-i) = f(t(i),u(:,i));
    end

    % Time stepping.
    for i = k:n
        ff = [f(t(i),u(:,i)), ff(:,1:k-1)];    % new value of du/dt
        u(:,i+1) = u(:,i) + h*(ff*sigma);      % advance one step
    end

    u = u.';   % conform to MATLAB output convention
end