function [ab4T, ab4W] = ab_4step(Fun, t0, tlast, w0, N)
%   function [ab4T, ab4W] = ab_4step(Fun, t0, tlast, w0, N)
%	Adams-Bashford four step method to solve
%	    y' = Fun(t,y), on [t0, tlast] with y(t0)=w0 with N+1 points
%
% INPUT:
% Fun   - String containing name of function m-file (enclose in '  ')
% t0    - Initial value of t.
% tlast - Last value of t.
% w0    - Initial value.
% N     - The number of steps
%
% OUTPUT:
% ab4T  - time points (column-vector).
% ab4W  - Solution column-vector per ab4T-values.
%
% The result can be displayed by: plot(ab4T, ab4W).

if nargin < 5, 
 error('Not enough arguments, enter help ab_4step for help')
end;
% Initialization
ab4T = t0;
h = (tlast - t0)/N;
ab4W = w0;
%% RK of order 4 to find first few steps
for i=1:3;
 k1= h*feval(Fun,ab4T(i),ab4W(i));  
 k2= h*feval(Fun,ab4T(i)+(h/2),ab4W(i)+(k1/2));  
 k3= h*feval(Fun,ab4T(i)+(h/2),ab4W(i)+(k2/2));  
 k4= h*feval(Fun,ab4T(i)+h,ab4W(i)+k3);  
 w = ab4W(i) + (k1+2*k2+2*k3+k4)/6; 
 t = ab4T(i) + h; clear k1 k2 k3 k4;
 ab4T=[ab4T;t]; ab4W = [ab4W;w];
end;

% The main loop
   for i=4:N;
      k_i = feval(Fun,ab4T(i),ab4W(i));
      k_i1 = feval(Fun,ab4T(i-1),ab4W(i-1));
      k_i2 = feval(Fun,ab4T(i-2),ab4W(i-2));
      k_i3 = feval(Fun,ab4T(i-3),ab4W(i-3));
      w=w+(h/24)*(55*k_i-59*k_i1+37*k_i2-9*k_i3);
      t=t+h;
      ab4T=[ab4T;t]; ab4W = [ab4W;w];
   end;

