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

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

% The main loop
   for i=2:N;
      k_i = feval(Fun,ab2T(i),ab2W(i));
      k_i1 = feval(Fun,ab2T(i-1),ab2W(i-1));
      w=w+(h/2)*(3*k_i-k_i1);
      t=t+h;
      ab2T=[ab2T;t]; ab2W = [ab2W;w];
   end;

