% the program takes the initial data y_0 at the initial time x_0 and 
% computes up to time x_final using N timesteps of size dx.  It creates
% two vectors, x and y, which you can then plot against each other
% plot(x,y)
%
% To call the routine, use the command:
%
% [x,y] = adams_bashforth_4(y_0,x_0,x_final,N,@f)
%
function [x,y] = adams_bashforth_4(y_0,x_0,x_final,N,f)

% step-size:
dx = (x_final-x_0)/N;

% for dy/dx = f(x,y(x))
%
% y(x+dx) ~ y(x)+dx/12*(23*f(x,y(x))-16*f(x-dx,y(x-dx))+5*f(x-2*dx,y(x-2*dx)))

% pux initial data into vector
x(1) = x_0;
y(:,1) = y_0;

% I need to get y(x_0+dx) with local truncation error O(h^4).  I'll
% do this with Richardson extrapolation on Euler timestepping...
ii=1;
F1 = feval(f,x(ii),y(:,ii));
F_now = F1;
% take one big step to reach y_big
y_big = y(:,ii) + dx*F1;
% take two small steps to reach y_med
y_temp = y(:,ii) + dx/2*F1;
F = feval(f,x(ii)+dx/2,y_temp);
y_med = y_temp + dx/2*F;
% take four smaller yet steps to reach y_small
% y at x + dx/4
y_temp = y(:,ii) + dx/4*F1;
F = feval(f,x(ii)+dx/4,y_temp);
% y at x + 2*(dx/4)
y_temp = y_temp + dx/4*F;
F = feval(f,x(ii)+2*(dx/4),y_temp);
% y at x + 3*(dx/4)
y_temp = y_temp + dx/4*F;
F = feval(f,x(ii)+3*(dx/4),y_temp);
y_small = y_temp + dx/4*F;
% use Richardson extrapolation
rich_big = 2*y_med - y_big;
rich_small = 2*y_small - y_med;
y(:,ii+1) = (4/3)*rich_small - (1/3)*rich_big;
x(ii+1) = x_0 + ii*dx;

ii=2;
F_old = F_now;
F1 = feval(f,x(ii),y(:,ii));
F_now = F1;
% take one big step to reach y_big
y_big = y(:,ii) + dx*F1;
% take two small steps to reach y_med
y_temp = y(:,ii) + dx/2*F1;
F = feval(f,x(ii)+dx/2,y_temp);
y_med = y_temp + dx/2*F;
% take four smaller yet steps to reach y_small
% y at x + dx/4
y_temp = y(:,ii) + dx/4*F1;
F = feval(f,x(ii)+dx/4,y_temp);
% y at x + 2*(dx/4)
y_temp = y_temp + dx/4*F;
F = feval(f,x(ii)+2*(dx/4),y_temp);
% y at x + 3*(dx/4)
y_temp = y_temp + dx/4*F;
F = feval(f,x(ii)+3*(dx/4),y_temp);
y_small = y_temp + dx/4*F;
% use Richardson extrapolation
rich_big = 2*y_med - y_big;
rich_small = 2*y_small - y_med;
y(:,ii+1) = (4/3)*rich_small - (1/3)*rich_big;
x(ii+1) = x_0 + ii*dx;

ii=3;
F_older = F_old;
F_old = F_now;
F1 = feval(f,x(ii),y(:,ii));
F_now = F1;
% take one big step to reach y_big
y_big = y(:,ii) + dx*F1;
% take two small steps to reach y_med
y_temp = y(:,ii) + dx/2*F1;
F = feval(f,x(ii)+dx/2,y_temp);
y_med = y_temp + dx/2*F;
% take four smaller yet steps to reach y_small
% y at x + dx/4
y_temp = y(:,ii) + dx/4*F1;
F = feval(f,x(ii)+dx/4,y_temp);
% y at x + 2*(dx/4)
y_temp = y_temp + dx/4*F;
F = feval(f,x(ii)+2*(dx/4),y_temp);
% y at x + 3*(dx/4)
y_temp = y_temp + dx/4*F;
F = feval(f,x(ii)+3*(dx/4),y_temp);
y_small = y_temp + dx/4*F;
% use Richardson extrapolation
rich_big = 2*y_med - y_big;
rich_small = 2*y_small - y_med;
y(:,ii+1) = (4/3)*rich_small - (1/3)*rich_big;
x(ii+1) = x_0 + ii*dx;

for ii=4:N
  F_oldest = F_older;
  F_older = F_old;
  F_old = F_now;
  F_now = feval(f,x(ii),y(:,ii));
  y(:,ii+1) = y(:,ii) + dx/24*(55*F_now - 59*F_old + 37*F_older - 9*F_oldest);
  x(ii+1) = x_0 + ii*dx;
end