function [extT,extW,extH]=extra56(Fun,t0,tlast,w0,hmax,hmin,TOL);
%[extT,extW,extH]=extra56(Fun,t0,tlast,w0,hmax,hmin,TOL); 
%   Extrapolation Algorithm 5.6 from Burden and Faires
%   to approximate the solution to the problem
%	    y' = Fun(t,y), on [t0, tlast] with y(t0)=w0 
% INPUT:
% Fun   - String containing name of function m-file (enclose in '  ')
% t0    - Initial value of t.
% tlast - Last value of t.
% w0    - Initial value.
% hmax  - maximum step size;  defaults to (tlast-t0)/10
% hmin  - minimum step size;  defaults to (tlast-t0)/200
% TOL   - local truncation tolerence; defaults to 1.e-6
%
% OUTPUT:
% extT  - time points (column-vector).
% extW  - Solution column-vector per extT-values.
% extH  - Column of step sizes h
% The result can be displayed by: plot(extT, extH).

if nargin < 4, 
 error('Not enough arguments, enter help extra56 for guidelines')
end;
if nargin < 5, 
 hmax=(tlast-t0)/10;
end;
if nargin < 7, 
 TOL = 1.e-6;
end;
if nargin < 6, 
 hmin=(tlast-t0)/200;
end;
% Initialization %% Steps 1 and  2
extT = t0;
h = hmax;
extW = w0;
extH=0;
FLAG=1;
y(1)=w0;
NK=[2 3 4 6 8 12 16 24];
%% Step3
Q=zeros(7,7);
  for i=1:7;
     for j=1:i;
       Q(i,j)=(NK(i+1)./NK(j)).^2;
     end;
  end;

%%Step 4
while FLAG==1
   k=1;NFLAG=0;  %STEP 5
   while (k<=8 & NFLAG==0)  %step6
       HK=h./NK(k);          %STEP 7
       t=t0;
       w2=w0;
       w3=w2+HK.*feval(Fun,t,w2);
       t=t0+HK;
       for j=1:NK(k)-1;  %%STEP8
            w1=w2;
            w2=w3;
            w3=w1+2*HK.*feval(Fun,t,w2);
            t=t0+(j+1).*HK;
       end;
       y(k)=(w3+w2+HK*feval(Fun,t,w3))/2;  %%STEP 9
       if k>=2;    %%STEP 10
            j=k;v=y(1); %STEP 11
            while j>=2   %% STEP 12
                  y(j-1)=y(j)+(y(j)-y(j-1))./(Q(k-1,j-1)-1);
                  j=j-1;
            end;
            if abs(y(1)-v)<=TOL; %%STEP 13
                 NFLAG=1;
            end;
       end;
       k=k+1;  %%STEP 14
   end;
   k=k-1;  %%STEP 15
   if NFLAG==0   %%STEP 16
       h=h/2;  %%STEP 17
       if h< hmin;  %%STEP 18
           message='hmin violated'
           FLAG=0;
       end;
   else
       w0=y(1);   %%STEP 19
       t0=t0+h;
       extT=[extT;t0];extW=[extW;w0];extH=[extH;h];
       if abs(t0-tlast)<.5*hmin;  %%STEP 20
             FLAG=0;
       elseif t0+h>tlast;
             h=tlast-t0;
       elseif (k<=3 & h < .5*hmax)
             h=2*h;
       end;
   end;
end;
           
