function [x,y]=num6(fcn,x0,xe,y0,y1,n);
%
% effectively 5 stages, 6th order, P-stable method 
%
% y' = f(y)
%
% y(x0) = y0
% y(x0+h) = y1
% n = # of steps


% The coefficients
a=[0,0,0,0,0,0;
   0,0,0,0,0,0;
   443/33540,1091/2250,443/33540,0,153664/628875,153664/628875;
   -10614137757603/305435106787328,-19694681828047/200441788829184,726410229286345/47609697270474752,-50074796005/99249733632,-1/232,65/129;
   26305365/1027704832,21485939/73407488,-735790025/29803440128,21377845385/44705160192,2/33,-8/19;
   -11628195/511393792,-356459/7375872,1870595/71300096,-598971725/1140801536,-1/24,19/39];
b=[443/33540,1091/2250,443/33540,0,153664/628875,153664/628875];
c=[-1,0,1,-21/37,15/28,-15/28]';


s=6;                         % stages
h=(xe-x0)/n;                 % step length
d=1e-8;                      % difference for evaluation of Jacobians
m=length(y0);                % dimension of system
x=[x0 x0+h zeros(1,n-1)]';   % output of x
y=zeros(m,n+1);              % output of y
id=eye(m);
aid=kron(a(3:s,:),id);
ids=eye(s*m);
y(:,1)=y0;
y(:,2)=y1;
tol=1e-11;                   % tolerance for Newton iteration
cc1=1/6*h^2*c(3:s).*(c(3:s)+1).*(c(3:s)+2); % used for initialization of z
cc0=1/6*h^2*c(3:s).*(c(3:s)+1).*(c(3:s)-1);
F=zeros(s*m,1);
f1=feval(fcn,y0);

for k=2:n,
    
   f0=f1;
   F(1:m)=f0;
   
   % forming the Jacobian
   jac=zeros(m);
   f1=feval(fcn,y(:,k));
   for o=1:m,
       jac(:,o)=(feval(fcn,y(:,k)+d*id(:,o))-f1)/d;
   end;
   F(m+1:2*m)=f1;
   aa=ids-h*h*kron(a,jac);
   part22=inv(aa(2*m+1:m*s,2*m+1:m*s));          % inversion of implicit part only

   % Modified Newton
   z=kron(cc1,f1)-kron(cc0,f0);  % Initial value of z
   
   for o=3:s,
       F((o-1)*m+1:o*m)=feval(fcn,(1+c(o))*y(:,k)-c(o)*y(:,k-1)+z((o-3)*m+1:(o-2)*m));
   end;
   dife=part22*(-z+h*h*aid*F);
   
   y(:,k+1)=2*y(:,k)-y(:,k-1)+z(1:m)+dife(1:m);  % use 3rd stage for b
   x(k+1)=x(k)+h;
end;