function sol = bvp4c(ode, bc, solinit, options, varargin);
%BVP4C  Solve two-point boundary value problems for ODEs by collocation.   
%   SOL = BVP4C(ODEFUN,BCFUN,SOLINIT) integrates a system of ordinary
%   differential equations of the form y' = f(x,y) on the interval [a,b],
%   subject to general two-point boundary conditions of the form
%   bc(y(a),y(b)) = 0. ODEFUN is a function of two arguments: a scalar X
%   and a vector Y. ODEFUN(X,Y) must return a column vector representing
%   f(x,y). BCFUN is a function of two vector arguments. BCFUN(YA,YB) must
%   return a column vector representing bc(y(a),y(b)). SOLINIT is a structure
%   with fields named   
%       x -- ordered nodes of the initial mesh with 
%            SOLINIT.x(1) = a, SOLINIT.x(end) = b
%       y -- initial guess for the solution with SOLINIT.y(:,i)
%            a guess for y(x(i)), the solution at the node SOLINIT.x(i)       
%
%   BVP4C produces a solution that is continuous on [a,b] and has a
%   continuous first derivative there. The solution is evaluated at points
%   XINT using the output SOL of BVP4C and the function BVPVAL:
%   YINT = BVPVAL(SOL,XINT). The output SOL is a structure with 
%       SOL.x  -- mesh selected by BVP4C
%       SOL.y  -- approximation to y(x) at the mesh points of SOL.x
%       SOL.yp -- approximation to y'(x) at the mesh points of SOL.x
%
%   SOL = BVP4C(ODEFUN,BCFUN,SOLINIT,OPTIONS) solves as above with default
%   parameters replaced by values in OPTIONS, a structure created with the
%   BVPSET function. See BVPSET for details.  
%
%   SOL = BVP4C(ODEFUN,BCFUN,SOLINIT,OPTIONS,P1,P2...) passes constant, known
%   parameters P1, P2... to the functions ODEFUN and BCFUN, and to all 
%   functions specified in OPTIONS. Use OPTIONS = [] as a place holder if
%   no options are set.   
%   
%   Some boundary value problems involve a vector of unknown parameters p
%   that must be computed along with y(x):
%       y' = f(x,y,p)
%       0  = bc(y(a),y(b),p) 
%   For such problems the field SOLINIT.parameters is used to provide a guess
%   for the unknown parameters. On output the parameters found are returned
%   in the field SOL.parameters. The solution SOL of a problem with one set 
%   of parameter values can be used as SOLINIT for another set. Difficult BVPs 
%   may be solved by continuation: start with parameter values for which you can 
%   get a solution, and use it as a guess for the solution of a problem with 
%   parameters closer to the ones you want. Repeat until you solve the BVP 
%   for the parameters you want.
%
%   The function BVPINIT forms the guess structure in the most common 
%   situations:  SOLINIT = BVPINIT(X,YINIT) forms the guess for an initial mesh X
%   as described for SOLINIT.x and YINIT either a constant vector guess for the
%   solution or a function that evaluates the guess for the solution
%   at any point in [a,b]. If there are unknown parameters, 
%   SOLINIT = BVPINIT(X,YINIT,PARAMS) forms the guess with the vector PARAMS of 
%   guesses for the unknown parameters.  
%
%   Example
%         solinit = bvpinit([0 1 2 3 4],[1 0]);
%         sol = bvp4c(@twoode,@twobc,solinit);
%     solve a BVP on the interval [0 4] with differential equations and 
%     boundary conditions computed by functions twoode and twobc, respectively.
%     This example uses [0 1 2 3 4] as an initial mesh, and [1 0] as an initial 
%     approximation of the solution components at the mesh points.
%         xint = linspace(0,4);
%         yint = bvpval(sol,xint);
%     evaluate the solution at 100 equally spaced points in [0 4]. The first
%     component of the solution is then plotted with 
%         plot(xint,yint(1,:));
%   For more examples see TWOBVP, MAT4BVP and SHOCKBVP.   
%
%   See also BVPSET, BVPGET, BVPINIT, BVPVAL, @.
 
%   BVP4C is a finite difference code that implements the 3-stage Lobatto
%   IIIa formula. This is a collocation formula and the collocation 
%   polynomial provides a C1-continuous solution that is fourth order
%   accurate uniformly in [a,b]. Mesh selection and error control are based
%   on the residual of the continuous solution. Analytical condensation is
%   used when the system of algebraic equations is formed.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2000 The MathWorks, Inc. 
%   $Revision: 1.9 $  $Date: 2000/06/02 00:11:22 $

true = 1;
false = ~true;

% check input parameters
if nargin <3 
  error('Not enough input arguments')
end

if ~isstruct(solinit)
  error('The initial profile must be provided as a structure')
elseif ~isfield(solinit,'x')
  error(['The field ''x'' not present in ''' inputname(3) ''''])
elseif ~isfield(solinit,'y')
  error(['The field ''y'' not present in ''' inputname(3) ''''])
end

x = solinit.x(:)';     % row vector
y = solinit.y;   

if isempty(x) | (length(x)<2)
  error(['''' inputname(3) '.x'' must contain at least the two end points'])
else
  N = length(x);      % number of mesh points
end
xdir = sign(x(end)-x(1));
if any( xdir * diff(x) <= 0 )
  error (['The entries in ''' inputname(3) ...
          '.x'' must strictly increase or decrease']);
end

if isempty(y)
  error(['No initial guess provided in ''' inputname(3) '.y'''])
end
if length(y(1,:)) ~= N
  error(['''' inputname(3) '.y'' not consistent with ''' ...
         inputname(3) '.x'''])
end

n = size(y,1);	% size of the DE system
nN = n*N;

% stats
nODEeval = 0;
nBCeval = 0;

% set the options
if nargin<4
  options = [];
end

% parameters
knownPar = varargin;
unknownPar = isfield(solinit,'parameters');
if unknownPar    
  vect_par = solinit.parameters(:);   % copy input to a column  vector
  npar = length(vect_par);  
  ExtraArgs = [{vect_par} knownPar];
else
  npar = 0;      
  vect_par = [];  
  ExtraArgs = knownPar;    
end               
  
% Get options and set the defaults
rtol = bvpget(options,'RelTol',1e-3);      
if (length(rtol) ~= 1) | (rtol<=0)
  error('RelTol must be a positive scalar.');  
end
if rtol < 100*eps
  rtol = 100*eps;
  warning(['RelTol has been increased to ' num2str(rtol) '.']);
end  
atol = bvpget(options,'AbsTol',1e-6);
if length(atol) == 1
  atol = ones(n,1)*atol;
else
  atol = atol(:);
end
if length(atol) ~= n
  error(sprintf(['Solving %s requires a scalar AbsTol, '...
                 'or a vector AbsTol of length %d'],funstring(ode),n));
end
if any(atol<=0)
  error('AbsTol must be positive')
end  

threshold = atol/rtol;

% analytical Jacobians
Fjac = bvpget(options,'FJacobian');
BCjac = bvpget(options,'BCJacobian');
Nmax = bvpget(options,'Nmax',floor(1000/n));  
printstats = strcmp(bvpget(options,'Stats','off'),'on');

yp = zeros(size(y));
for i=1:N
  yp(:,i) = feval(ode,x(i),y(:,i),ExtraArgs{:});  
end  
nODEeval = nODEeval + N;

maxNewtIter = 4;  
maxProbes = 4;    % weak line search 
needGlobJac = true;

done = false;
     
% THE MAIN LOOP:
while ~done       
  Y = zeros(nN+npar,1);
  Y(1:nN)  = reshape(y, nN,1);      
  if unknownPar
    Y(nN+1:nN+npar) = vect_par;
  end  
  
  [RHS,Fmid,NF] = colloc_eqns(n,x,Y,yp,ode,bc,npar,ExtraArgs); 	
  nODEeval = nODEeval + NF;
  nBCeval = nBCeval + 1;
    
  for iter=1:maxNewtIter	
    if needGlobJac
      % setup and factor the global Jacobian            
      [dPHIdy,NF,NBC] = colloc_Jac(n,x,Y,yp,ode,bc,Fjac,BCjac,npar,ExtraArgs); 
      needGlobJac = false;                     
      nODEeval = nODEeval + NF;                
      nBCeval = nBCeval + NBC;      

      % explicit row scaling              
      RowScale = 1 ./ max(abs(dPHIdy)')';
      sizeRSc = length(RowScale);
      scalMatrix = spdiags(RowScale,0,sizeRSc,sizeRSc);
      dPHIdy = scalMatrix  * dPHIdy;      
      % check the Jacobian for singularity
      singJac = (condest(dPHIdy) > 1e10);
      if singJac           
        warning('an ill-conditioned Jacobian -- the mesh will be halved'); 
        break     
      end                     
      [L,U,P] = lu(dPHIdy);
      scalMatrix = P * scalMatrix;                       
    end
                    
    % find the Newton direction, 
    delY = U\(L\( scalMatrix * RHS ));    
    distY = norm(delY);        

    % weak line search with an affine stopping criterion
    lambda = 1;
    for probe=1:maxProbes     
      Ynew = Y - lambda*delY;                

      y  = reshape(Ynew(1:nN),n,N);
      if unknownPar
        vect_par = Ynew(nN+1:nN+npar);   
        ExtraArgs = [{vect_par} knownPar];
      end
      yp = zeros(n,N);
      for i=1:N
        yp(:,i) = feval(ode,x(i),y(:,i),ExtraArgs{:});  
      end  
      nODEeval = nODEeval + N; 
      
      [RHS,Fmid,NF] = colloc_eqns(n,x,Ynew,yp,ode,bc,npar,ExtraArgs);
      nODEeval = nODEeval + NF;
      nBCeval = nBCeval + 1;

      distYnew = norm(U \ (L \ (scalMatrix * RHS)));

      if (distYnew < 0.9*distY)       
        break		
      else
        lambda = 0.5*lambda;		
      end
    end  
    
    Y  = Ynew;   % y, yp, vect_par, ExtraArgs, and RHS 
                 % are consistent with Ynew 
        
    needGlobJac = (distYnew > 0.1*distY);
          
    if distYnew < 0.1*rtol
      break
    end        
  end  
  
  if singJac       
    NN = 2*N-1;   % halve the mesh 
    if NN > Nmax
      msg = sprintf(...
        [ 'Unable to refine the mesh any further -- '...
          'the Jacobian of the collocation equations is singular']);
      error(msg)    
    end    
    xx = zeros(1,NN);
    xx(1:2:NN) = x;
    xx(2:2:NN-1) = x(1:N-1) + 0.5*diff(x);
    
    % a new solution profile 
    yy = zeros(n,NN);
    yy(:,1:2:NN) = y;
    
    temp = struct('x',x,'y',y,'yp',yp);
    yy(:,2:2:NN-1) = bvpval(temp,xx(2:2:NN-1));
        
    yyp = zeros(size(yy));
    yyp(:,1:2:NN) = yp;
    for i=2:2:NN-1 
      yyp(:,i) = feval(ode,xx(i),yy(:,i),ExtraArgs{:}); 
    end
    nODEeval = nODEeval + (NN-1)/2;
    N = NN;
    nN = n*N;  
    x = xx;
    y = yy;
    yp = yyp;
    needGlobJac = true;    
    
  else
   
    % residual at the midpoints is related to the error
    % in satisfying the collocation equations
    h = diff(x);            
    NewtRes = reshape(RHS(n+npar+1:end),n,N-1); 

    for i=1 : N-1
      temp = 1.5 * NewtRes(:,i) / h(i);   
      res_mid(i) = norm (temp ./ max(abs(Fmid(:,i)),threshold));
    end
    [res,NF] = residual(ode,x,y,yp,res_mid,threshold,ExtraArgs);
    nODEeval = nODEeval + NF;     
  
    maxres = max(res);
    if maxres < rtol      
      done = true;
    else 	
      % redistribute the mesh
      [N,x,y] = new_profile(n,x,y,yp,res,rtol,Nmax);
      if N>Nmax    
        msg = sprintf(...
          [ 'Unable to meet the tolerance without using more than %d '...
            'mesh points. \n The last mesh of %d points and ' ...
            'the solution are available in the output argument. \n ',...
            'The maximum residual is %e, while requested accuracy %e \n'],...
          Nmax,length(x),max(res),rtol);       
        warning(msg); 
        sol.x = x;
        sol.y = y;
        sol.yp = yp;
        if unknownPar
          sol.parameters = vect_par;
        end          
        return  
      end      
      nN = n*N;           
      yp = zeros(n,N);
      for i=1:N
        yp(:,i) = feval(ode,x(i),y(:,i),ExtraArgs{:});    
      end
      nODEeval = nODEeval + N;    
      needGlobJac = true;	
    end  
  end     
end    % while

% Output
sol.x = x;
sol.y = y;
sol.yp = yp;
if unknownPar
  sol.parameters = vect_par;
end  

% Stats
if printstats
  fprintf('The solution was obtained on a mesh of %g points.\n',N);
  fprintf('The maximum residual is %10.3e. \n',maxres); 
  fprintf('There were %g calls to the ODE function. \n',nODEeval); 
  fprintf('There were %g calls to the BC function. \n',nBCeval); 
end

%--------------------------------------------------------------------------

function [Sx,Spx] = Hermite_cubic (xx, xn,yn,ypn,xnp1,ynp1,ypnp1);
%HERMITE_CUBIC Evaluate the cubic Hermite interpolant and its first
%  derivative 

h = xnp1 - xn;
t = xx - xn;
slope = (ynp1-yn)/h;
c = (3*slope - 2*ypn - ypnp1)/h;
d = (ypn+ypnp1-2*slope) / h^2;
Sx = ((d*t+c)*t+ypn)*t+yn;
if nargout>1
  Spx = (3*d*t+2*c)*t+ypn;
end

%---------------------------------------------------------------------------

function [res,nfcn] = residual (Fcn, x, y, yp, mid_res, threshold, ExtraArgs)  
%RESIDUAL  Compute L2-norm of the residual using 5-point Lobatto quadrature

N = length(x);   
res = zeros(1,N-1);    
res_sample = zeros(3,1);   
% Lobatto quadrature
lob4 = sqrt(3/7);  
lob2 = -lob4;                               
lobw24 = 49/90;   
lobw3 = 32/45;    
 
h = diff(x);   % interval lengths

for i=1:N-1  
  h2 = h(i)/2; 
  xmid = (x(i)+x(i+1))/2; 
  
  % discrete norm of the residual at Lob2 
  x2 = xmid + lob2*h2;
  [y2,yp2] = Hermite_cubic(x2,x(i),y(:,i),yp(:,i),x(i+1),y(:,i+1),yp(:,i+1));
  f2 = feval(Fcn,x2,y2,ExtraArgs{:});
  res_sample(1) = norm( (yp2-f2) ./ max(abs(f2),threshold) );  

  % discrete norm of the residual at Lob3, the midpoint
  res_sample(2) = mid_res(i);
  
  % discrete norm of the residual at Lob4    
  x4 = xmid + lob4*h2;
  [y4,yp4] = Hermite_cubic(x4,x(i),y(:,i),yp(:,i),x(i+1),y(:,i+1),yp(:,i+1));
  f4 = feval(Fcn,x4,y4,ExtraArgs{:});
  res_sample(3) = norm( (yp4-f4) ./ max(abs(f4),threshold) );
    
  % the residual on [x(i),x(i+1)]
  temp = abs(h2) * (lobw24*(res_sample(1)^2 + res_sample(3)^2) +...
                    lobw3*res_sample(2)^2);
  res(i) = sqrt(temp);
end
nfcn = 2*(N-1);

%---------------------------------------------------------------------------

function [NN,xx,yy] = new_profile(n,x,y,yp,res,rtol,Nmax)
%NEW_PROFILE  Redistribute mesh points and approximate the solution

N = length(x);
last_int = N-1;
h = diff(x);
xx(1) = x(1);
yy(:,1) = y(:,1);
NN = 1;
i=1;
removed = 0;
while i<=last_int
  if res(i) > rtol     % introduce points   
    if res(i) > 100*rtol   
      Ni = 2;           
    else
      Ni = 1;     
    end  
    hi = h(i) / (Ni+1);
    for j=1:Ni
      xx(NN+j) = xx(NN+j-1) + hi;        
      yy(:,NN+j) = Hermite_cubic(xx(NN+j),x(i),y(:,i),yp(:,i),...
                           x(i+1),y(:,i+1),yp(:,i+1));
    end
    NN = NN + Ni;
  else    
    if i <= last_int-2    
      if max(res(i+1:i+2)) < rtol   % try to remove points        
        hnew = (h(i)+h(i+1)+h(i+2))/2;
        C1 = res(i)/(h(i)/hnew)^(7/2);
        C2 = res(i+1)/(h(i+1)/hnew)^(7/2);
        C3 = res(i+2)/(h(i+2)/hnew)^(7/2);
        pred_res = max([C1,C2,C3]); 

        if pred_res < 0.5 * rtol   % replace 3 intervals with 2 
          removed = removed+1;
          xx(NN+1) = xx(NN) + hnew;
          yy(:,NN+1) = Hermite_cubic(xx(NN+1),x(i+1),y(:,i+1),yp(:,i+1),...
                               x(i+2),y(:,i+2),yp(:,i+2));
          NN = NN+1;  
          i = i+2;    
        end        
      end  
    end      
  end
  NN = NN+1;
  xx(NN) = x(i+1);   % preserve the next mesh point
  yy(:,NN) = y(:,i+1);
  i = i+1;
end    
if NN > Nmax  % too many mesh points requested
  xx = x; 
  yy = y;     % return the previous solution 
  return
end    

%---------------------------------------------------------------------------

function [Phi,Fmid,nfcn] = colloc_eqns(n, x, Y, F, Fcn, Gbc, npar, ExtraArgs)
%COLLOC_EQNS  Evaluate the system of collocation equations Phi(Y).  
%   The derivative approximated at the midpoints and returned in Fmid is
%   used to estimate the residual. 

N = length(x);
nN = n*N;              
Phi = zeros(nN+npar,1);
Fmid = zeros(n,N-1);    

h = diff(x);
xip05 = x(1:N-1) + 0.5*h;

% Boundary conditions
yip1 = Y(1:n);         % y(a)
yN = Y(nN-n+1:nN);     % y(b)
Phi(1:n+npar) = feval(Gbc,yip1,yN,ExtraArgs{:});    

Fip1 = F(:,1);

rowsY = (n+1:2*n);

for i = 1:N-1
  yi = yip1;
  Fi = Fip1;
  yip1 = Y(rowsY);
  Fip1 = F(:,i+1);
  
  yip05 = (yi+yip1)/2-h(i)/8*(Fip1-Fi);  
  Fip05 = feval(Fcn,xip05(i),yip05,ExtraArgs{:});  % exact derivative 

  % the Lobatto IIIa formula
  Phi(rowsY+npar) = (yip1 - yi - h(i)/6 * (Fi + 4*Fip05 + Fip1));

  % the derivative at the midpoint
  Fmid(:,i) = Fip05;    % will be used to estimate the residual

  rowsY = rowsY+n;
end
nfcn = N-1;

%---------------------------------------------------------------------------

function [dPHIdy,nfcn,nbc] = colloc_Jac(n, x, Y, F, ode, bc, ...
                                          Fjac, BCjac, npar, ExtraArgs); 
%colloc_Jac  Form the Jacobian of the system of collocation eqns

nfcn = 0;   
nbc = 0;
N = length(x);
nN = n*N;               
In = eye(n);

ya = Y(1:n);          
yb = Y(nN-n+1:nN);   
                               
rows = (n+1:2*n);   % define a region of 
cols = (1:n);       % the global Jacobian

if npar == 0    % no unknown parameters

  dPHIdy = spalloc(nN,nN,2*nN*n);      % storage
  
  % BC
  if isempty(BCjac)   % use numerical approx
    [dGdya,dGdyb,nbc] = BCnumjac(bc,ya,yb,n,npar,ExtraArgs);        
  else  % use analytical Jacobian 
    [dGdya,dGdyb] = feval(BCjac,ya,yb,ExtraArgs{:});
  end   
  dPHIdy(1:n,1:n) = dGdya;  % and postpone storing dGdyb        
  
  % Collocation equations
  if isempty(Fjac)  % use numerical approx
    
    xip1 = x(1);
    yip1 = Y(cols);
    Fip1 = F(:,1);
    threshold = repmat(1e-6,n,1);
    FAC = [];
    [Jip1,FAC,nF] = Fnumjac(ode,xip1,yip1,Fip1,n,npar,threshold,FAC,ExtraArgs); 
    nfcn = nfcn+nF;
        
    for i = 1:N-1 
      % the current grid point
      xi = xip1;            
      yi = yip1;
      Fi = Fip1;
      Ji = Jip1;
      % the next grid point
      xip1 = x(i+1);
      yip1 = Y(cols+n);
      Fip1 = F(:,i+1);
      [Jip1,FAC,nF] = Fnumjac(ode,xip1,yip1,Fip1,n,npar,threshold,FAC,ExtraArgs);
      nfcn = nfcn+nF;
      % the midpoint
      hi = xip1-xi;
      xip05 = (xi+xip1)/2; 
      yip05 = (yi+yip1)/2-hi/8*(Fip1-Fi);  
      Fip05 = feval(ode,xip05,yip05,ExtraArgs{:});
      [Jip05,FAC,nF] = Fnumjac(ode,xip05,yip05,Fip05,n,npar, ...
                               threshold,FAC,ExtraArgs);
      nfcn = nfcn+nF+1;
      % assembly
      twiceJip05 = 2*Jip05;
      dPHIdy(rows,cols) = -(In+hi/6*(Ji+twiceJip05*(In+hi/4*Ji)));
      cols = cols + n;
      dPHIdy(rows,cols) = In-hi/6*(Jip1+twiceJip05*(In-hi/4*Jip1));
      rows = rows+n;   % next equation
    end        
  else   % use analytical Jacobian          
    xip1 = x(1);
    yip1 = Y(cols);
    Fip1 = F(:,1);
    Jip1 = feval(Fjac,xip1,yip1,ExtraArgs{:});
    
    for i = 1:N-1 
      % the current grid point
      xi = xip1;            
      yi = yip1; 
      Fi = Fip1;
      Ji = Jip1;
      % the next grid point
      xip1 = x(i+1);
      yip1 = Y(cols+n);
      Fip1 = F(:,i+1);
      Jip1 = feval(Fjac,xip1,yip1,ExtraArgs{:});    
      % the midpoint
      hi = xip1-xi;
      xip05 = (xi+xip1)/2; 
      yip05 = (yi+yip1)/2-hi/8*(Fip1-Fi);  
      Jip05 = feval(Fjac,xip05,yip05,ExtraArgs{:});  % recompute the Jacobian
      % assembly
      twiceJip05 = 2*Jip05;
      dPHIdy(rows,cols) = -(In+hi/6*(Ji+twiceJip05*(In+hi/4*Ji)));
      cols = cols + n;
      dPHIdy(rows,cols) = In-hi/6*(Jip1+twiceJip05*(In-hi/4*Jip1));
      rows = rows+n;   % next equation
    end
  end
  dPHIdy(1:n,nN-n+1:nN) = dGdyb;    % computed earlier
  
else  % there are unknown parameters

  dPHIdy = spalloc(nN+npar,nN+npar,(nN+npar)*(2*n+npar));  % storage
  last_cols = zeros(nN+npar,npar);   % accumulator
  
  % BC
  if isempty(BCjac)   % use numerical approx
    [dGdya,dGdyb,nbc,dGdpar] = BCnumjac(bc,ya,yb,n,npar,ExtraArgs);        
  else  % use analytical Jacobians 
    [dGdya,dGdyb,dGdpar] = feval(BCjac,ya,yb,ExtraArgs{:});
  end   
  dPHIdy(1:n+npar,1:n) = dGdya;  % and postpone storing dGdyb        
  last_cols(1:n+npar,:) = dGdpar;
  
  % accomodate additional npar of boundary conditions
  rows = rows+npar;

  if isempty(Fjac)  % use numerical approx    
    xip1 = x(1);
    yip1 = Y(cols);
    Fip1 = F(:,1);    
    threshold = repmat(1e-6,n+npar,1);
    FAC = [];
    [Jip1,FAC,nF,dFdpar_ip1] = Fnumjac(ode,xip1,yip1,Fip1,n, ...
                                       npar,threshold,FAC,ExtraArgs); 
    nfcn = nfcn+nF;
    
    for i = 1:N-1
      % the current grid point   
      xi = xip1;            
      yi = yip1; 
      Fi = Fip1;
      Ji = Jip1;
      dFdpar_i = dFdpar_ip1;     
      % the next grid point
      xip1 = x(i+1);
      yip1 = Y(cols+n);
      Fip1 = F(:,i+1);
      [Jip1,FAC,nF,dFdpar_ip1] = Fnumjac(ode,xip1,yip1,Fip1,n,npar,...
                                         threshold,FAC,ExtraArgs);
      nfcn = nfcn+nF;      
      % the midpoint 
      hi = xip1-xi;
      xip05 = (xi+xip1)/2; 
      yip05 = (yi+yip1)/2-hi/8*(Fip1-Fi);  
      Fip05 = feval(ode,xip05,yip05,ExtraArgs{:});
      [Jip05,FAC,nF,dFdpar_ip05] = Fnumjac(ode,xip05,yip05,Fip05,n,...
                                           npar,threshold,FAC,ExtraArgs);
      nfcn = nfcn+nF+1;
      % assembly
      twiceJip05 = 2*Jip05;
      dPHIdy(rows,cols) = -(In+hi/6*(Ji+twiceJip05*(In+hi/4*Ji)));
      cols = cols + n;
      dPHIdy(rows,cols) = In-hi/6*(Jip1+twiceJip05*(In-hi/4*Jip1));
      last_cols(rows,:) = -hi*dFdpar_ip05 + hi^2/12*Jip05*...
                                            (dFdpar_ip1-dFdpar_i);
      rows = rows+n;   % next equation
    end
  else   % use analytical Jacobians      
    xip1 = x(1);
    yip1 = Y(cols);
    Fip1 = F(:,1);
    [Jip1,dFdpar_ip1] = feval(Fjac,xip1,yip1,ExtraArgs{:});    
    for i = 1:N-1
      % the current grid point   
      xi = xip1;            
      yi = yip1; 
      Fi = Fip1;
      Ji = Jip1;
      dFdpar_i = dFdpar_ip1;
      % the next grid point
      xip1 = x(i+1);
      yip1 = Y(cols+n);
      Fip1 = F(:,i+1);
      [Jip1, dFdpar_ip1] = feval(Fjac,xip1,yip1,ExtraArgs{:});    
      % the midpoint 
      hi = xip1-xi;
      xip05 = (xi+xip1)/2; 
      yip05 = (yi+yip1)/2-hi/8*(Fip1-Fi);  
      [Jip05, dFdpar_ip05] = feval(Fjac,xip05,yip05,ExtraArgs{:});        
      % assembly
      twiceJip05 = 2*Jip05;
      dPHIdy(rows,cols) = -(In+hi/6*(Ji+twiceJip05*(In+hi/4*Ji)));
      cols = cols+n;
      dPHIdy(rows,cols) = In-hi/6*(Jip1+twiceJip05*(In-hi/4*Jip1));
      last_cols(rows,:) = -hi*dFdpar_ip05 + hi^2/12*Jip05* ...
                                            (dFdpar_ip1-dFdpar_i); 
      rows = rows+n;   % next equation
    end    
  end
  dPHIdy(1:n+npar,nN-n+1:nN) = dGdyb;    % computed earlier
  dPHIdy(:,nN+1:nN+npar) = last_cols;  % accumulated
        
end

ndfdy = 2*N-1;

%---------------------------------------------------------------------------

function [dBCdya,dBCdyb,nbc,dBCdpar] = BCnumjac(bc,ya,yb,n,npar,ExtraArgs)
if npar > 0   % unknown parameters    
  argvect = [ya; yb; ExtraArgs{1}];   % stack column vectors 
else
  argvect = [ya; yb];        % stack column vectors
end
BCval = feval(bc,ya,yb,ExtraArgs{:});
thresh = repmat(1e-6,size(argvect));
fac = [];
[temp,fac,dummy,nbc] = numjac(@BCaux,bc,argvect,BCval,...
                              thresh,fac,0,[],[],n,npar,ExtraArgs);
dBCdya = temp(:,1:n);
dBCdyb = temp(:,n+1:2*n);
if npar>0  %unknown parameters
  dBCdpar = temp(:,2*n+1:end);
end
nbc = nbc+1;
  
%---------------------------------------------------------------------------

function BCval = BCaux(bc,argvec,n,npar,ExtraArgs);
ya = argvec(1:n);
yb = argvec(n+1:2*n);
if npar > 0
  ExtraArgs{1} = argvec(2*n+1:end);  % unknown parameters
end  
BCval = feval(bc,ya,yb,ExtraArgs{:});
  
%---------------------------------------------------------------------------

function [Ji,FAC,nF,dFdpar_i] = Fnumjac(ode,xi,yi,Fi,n,npar,thresh,FAC,ExtraArgs);
% compute the Jacobian of Faux(ode,argvect,xi,n,npar) wrt argvect
if npar >0
  argvect = [yi; ExtraArgs{1}];   % unknown parameters - stack column vectors  
else
  argvect = yi;   
end    
  
[temp,FAC,dummy,nF] = numjac(@Faux,ode,argvect,Fi,...
                             thresh,FAC,0,[],[],xi,n,npar,ExtraArgs);  
Ji = temp(:,1:n);
if npar>0  %unknown parameters
  dFdpar_i = temp(:,n+1:end);
end
  
%---------------------------------------------------------------------------

function Fval = Faux(ode,argvect,x,n,npar,ExtraArgs);
y = argvect(1:n);     
if npar > 0
  ExtraArgs{1} = argvect(n+1:end);    % unknown parameters      
end
Fval = feval(ode,x,y,ExtraArgs{:});       


  
