%--------------------------------------------------------------------------------------------------------
function [xout,yout,yend,ireject]=num6v(fcn,x0,xe,y0,dy0,tol);
% Variable step Numerov 6th order (CMA2003)

f0=feval(fcn,x0,y0);                             % f0
h=(.5*tol/(max(abs(f0))+tol))^.25;               % initial step
[x1,y1]=rkn86(fcn,x0,x0+h,y0,dy0,tol/100);       % initial y1
y1=y1(end,:)';
x=x0+h;                                          % advance x
f1=feval(fcn,x,y1);                              % f1
xout = [x0 x0+h zeros(1,634998)];                % preallocate output for x
m=length(y0);                                    % dimension of system
yout=zeros(m,635000);                            % preallocate output of y
yout(:,1)=y0;
yout(:,2)=y1;
k=2;ireject=0;

% The main loop
   while (x < xe),

      ya1=(0.06366100187501753*f0+0.43633899812498247*f1)*h^2-0.61803398874989485*y0 ...
                                                                +1.61803398874989485*y1;
      fa1=feval(fcn,0.61803398874989485*h+ x,ya1);              % 1st stage

      ya2=(0.43633899812498247*f0+0.06366100187501753*f1)*h^2+1.61803398874989485*y0 ...
                                                                -0.61803398874989485*y1;
      fa2=feval(fcn,-1.61803398874989485*h+x,ya2);              % 2nd stage

      ya3=0.5*y0+0.5*y1+ ...
          (-0.03645833333333333*f0-0.114583333333333333*f1+0.03049011440755044*fa1   ...
                                                         -0.004448447740883774*fa2)*h^2;
      fa3=feval(fcn,x-1/2*h,ya3);                               % 3rd stage

      y2=-y0+2*y1+ ...
         (0.03333333333333333*f0+0.533333333333333333*f1+0.2218033988749895*fa1      ...
          -0.001803398874989485*fa2+0.213333333333333333*fa3)*h^2;
      f2=feval(fcn,x+h,y2);                                     % Compute 4th stage (FSAL)

                                                                % error estimation
      delta=max(abs(h^2*(-67*f0/480-67*f1/480+67*fa1/2400+67*fa2/2400+67*fa3/300)));

      % Update the solution only if the error is acceptable
      if delta <= tol*8,
         x = x + h;
         k=k+1;
         xout(k)=x;
         yout(:,k)=y2;
         if delta<tol/8,
             h=h+h;                              % Double the step
         else
             y0=y1;f0=f1;                        % Keep the same step
         end;
         y1=y2;
         f1=f2;
      else                                       % Halve the step
         ireject=ireject+1;
         fa4=feval(fcn,x-h/2,1/2*y0+1/2*y1+h^2*(-23/1920*f0-23/1920*f1+1/3200*fa1    ...
                                                               +1/3200*fa2-61/600*fa3));
         y0=1/2*y1+1/2*y0+h^2*(-23/1920*f0-23/1920*f1+1/3200*fa1+1/3200*fa2-61/600*fa4);
         f0=fa4;
         h=h/2;
      end;
   end;
   xout=xout(1:k);yout=yout(:,1:k);

   % since the end-point xe is not hitted, return yend by interpolation
   for j2=1:m,
       yend(j2)=polyval(polyfit(10*(xout(end-6:end)-xe),yout(j2,end-6:end),6),0);
   end;
return
%--------------------------------------------------------------------------------------------------------