function [t,u] = ab4(dudt,tspan,u0,n)
%AB4     4th-order Adams-Bashforth formula for an IVP.
% Input:
%   dudt    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(dudt,[a, a+(k-1)*h],u0,k-1);
u(:,1:k) = us(1:k,:).';

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

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

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