% -----------------------------------------------------------------------------------------------
function [x,y]=num4(fcn,x0,xe,y0,y1,n);
% 4th order, P-stable, Numerov type method with only 2 stages

h=(xe-x0)/n;           % the step-length
d=1e-6;                % difference for evaluation of Jacobians
tol=1e-4;              % tolerance for Newton
m=length(y0);          % dimension of system
x=(x0:h:xe)';          % output of x
y=zeros(m,n+1);        % output of y
y(:,1)=y0; 
y(:,2)=y1;
f0=feval(fcn,x(1),y0); 
f1=feval(fcn,x(2),y1);
irep=10000;

for k=2:n,

    if irep>3, % new Jacobian after Newton semi-failure
       jac=numjac(fcn,x(k),y1,f1,d*ones(m,1));
       [L1,U1,P1]=lu(eye(m)-0.110360360360360*jac*h^2+0.012008349813228*jac^2*h^4);
       ulp=U1\(L1\(P1));
    end;

    y2=-y0+2*y1+h^2*f1; f2=feval(fcn,x(k+1),y2);
    dife=10*tol; irep=0;

    while dife>tol,
          u3=(0.030020874533070*(f0+f2) + 0.007525818501428*f1)*h^2 + ...
              0.932432432432432*y0 + 0.135135135135135*y1 - 0.067567567567568*y2;
          f3=feval(fcn,x(k-1),u3);
          dife=ulp*((0.483333333333333*f0+ 0.833333333333333*f1-0.4*f3+0.0833333333333333*f2)*h^2 ...
                     -y0+2*y1-y2);
          y2=y2+dife; f2=feval(fcn,x(k+1),y2);
          irep=irep+1;
    end;
    y0=y1; f0=f1;
    y1=y2; f1=f2;
    y(:,k+1)=y1;
end;
% -----------------------------------------------------------------------------------------------