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

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

% The main loop
   for i=5:N;
      k_i = feval(Fun,ab5T(i),ab5W(i));
      k_i1 = feval(Fun,ab5T(i-1),ab5W(i-1));
      k_i2 = feval(Fun,ab5T(i-2),ab5W(i-2));
      k_i3 = feval(Fun,ab5T(i-3),ab5W(i-3));
      k_i4 = feval(Fun,ab5T(i-4),ab5W(i-4));
      w=w+(h*(1901*k_i-2774*k_i1+2616*k_i2-1274*k_i3+251*k_i4)/720);
      t=t+h;
      ab5T=[ab5T;t]; ab5W = [ab5W;w];
   end;

