function [ab3T, ab3W] = ab_3step(Fun, t0, tlast, w0, N)
%  function [ab3T, ab3W] = ab_3step(Fun, t0, tlast, w0, N)
%	Adams-Bashford three 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:
% ab3T  - time points (column-vector).
% ab3W  - Solution column-vector per ab3T-values.
%
% The result can be displayed by: plot(ab3T, ab3W).

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

% The main loop
   for i=3:N;
      k_i = feval(Fun,ab3T(i),ab3W(i));
      k_i1 = feval(Fun,ab3T(i-1),ab3W(i-1));
      k_i2 = feval(Fun,ab3T(i-2),ab3W(i-2));
      w=w+(h/12)*(23*k_i-16*k_i1+5*k_i2);
      t=t+h;
      ab3T=[ab3T;t]; ab3W = [ab3W;w];
   end;

