Problem with Newton's methode for ODE system

Hello! I have got a system of 3 ODEs in disret form and want to solve them with Newton's method. I found 2 solvers here in the fileexchange implementing this method but with both has got an error "Index exceeds array bounds.". Any help would be appreciated.

function Gas_system_4v5_3
clc
close all
n =100;
x0 = [1000*ones(n,1);1*ones(n,1);1*ones(n,1)];
x = ((1:n)-1)/(n-1);
 %  some parameters
cfm=5*10^2; 
    l=1;
Pe_gm=l*cfm/1e+4; 
 %  solvers choise
Solver=1;
switch Solver
   case 1
[sol, error] = newtonraphson(@(x)fun(x,n),x0, 1000) ;       
case 2
    
  %  options1 = ['tolfun',1e-5, 'tolx',1e-4 ,'maxiter' , 1000 ];
sol = newton(@(x)fun(x,n), [1000*ones(n,1);1*ones(n,1);1*ones(n,1)] ,'maxiter',100,'tolfun',1e-3);
end

[~,Pe_g, Pe_dk] = fun(sol,n);
 

norm(fun(sol,n))
figure
plot(x*l*(Pe_gm)^-0.5,sol(1:n))
hold on
figure
plot(x*l*(Pe_gm)^-0.5,sol(n+1:2*n))
hold on
figure
plot(x*l*(Pe_gm)^-0.5,sol(2*n+1:3*n))
hold on

    function [res, Pe_g, Pe_dk] = fun(Yg,n)
  
    x = ((1:n)-1)/(n-1);
    dx = 1/(n-1);
% initial conditions_______________________________________________________
    T_fout =750;
    Y_gHMXout =0.8;
    Y_HMXprod0=0.2;
%-vectors_________________________________________________________________
    T_g_ = Yg(1:n);
    Y_gHMX_ = Yg(n+1:2*n);
    Y_HMXprod_ = Yg(2*n+1:3*n);

    
% paramters________________________________________________________________

    Y_gTf=0;
    Y_BF3 =0;
    Y_CF2=0;
    Y_C=0;
    Y_C2=0;
    Y_cTf=0.2;
    Y_B=0.8;
    fg=0.3;
    Q_gHMX=1359.8;
    %
    P =20;
    
    m=2;
    M_Tf = 100.014;
    
    R=8.314;
    ro_cTf=2.2;
    M_N2O=44.013;
    M_NO2=46.005;
    M_NO=30.006;
    M_CH2O=30.0258;
    M_H2O=18.0148;
    M_CO=28.01;
    TmTf=603.15;
    kTf=8*10^14; % C. D.  Doyle, J. Appl. Polymer Sci., 5,285 (1961).
    ETf=280328;
    ro_B =2.34;
    A=1;
    %
    L_g=0.1;
    
    
    Y_NO2_f =0.1;
    Y_NO_f=0.2;
    Y_N2O_f=0.3;
    Y_H2O_f=0.1;
    Y_CO_f=0.1;
    Y_CH2O_f=0.1;
    
    M_B= 10.8;
    M_HMX=296.1552;
    M_C2F4=100;
    M_BF3= 67.8;
    dH_BF3 = -1135.62;
    M_CF2=50;
    dH_CF2=- 180; %кДж/моль
    M_C=12;
    dH_C=716.68 ;
    
    M_C2=24;
    dH_C2=837.74;
    D = 0.5; % см2/с
    wm=m/l;
    
    dHm=cfm;
% functions for ODE________________________________________________________
    tau =@(T_g)(T_g/1000);
    
    Pn2=@(k0,k1, T_g)(k0+k1*tau(T_g));
    Pn5=@(k0,k1,k2,k3,k4, T_g)(k0+k1*tau(T_g)+k2*tau(T_g).^2+k3*tau(T_g).^3+k4./tau(T_g).^2);
    Ar=@(k,E)(k*exp(-E./(R*T_g_)));
    
    M_HMXprod=M_NO2*Y_NO2_f+M_NO*Y_NO_f+M_N2O*Y_N2O_f+M_H2O*Y_H2O_f+M_CO*Y_CO_f+M_CH2O*Y_CH2O_f;
    Mg=M_HMXprod*Y_HMXprod_+M_HMX*Y_gHMX_+M_Tf*Y_gTf + M_BF3*Y_BF3 + M_CF2*Y_CF2 +M_C*Y_C;
    
    % плотности
    rog=Mg.*P./(R.*T_g_);
    roc = ro_cTf*Y_cTf + ro_B*Y_B ;
    ro_gf=(1-fg)*roc+fg*rog;
    
    % теплоемкости газов
    c_NO2(T_g_<1200,1)=Pn5(16.108,75.89,-54.3874,14.30777,0.239423,T_g_(T_g_<1200))/M_NO2;
    c_NO2(T_g_>=1200,1)=Pn5(56.82541,0.738053,-0.144721, 0.009777,-5.459911,T_g_(T_g_>=1200))/M_NO2;
    c_NO(T_g_<1200,1)=Pn5(23.83, 12.589, -1.139, -1.497459, 0.214194,T_g_(T_g_<1200))/M_NO;
    c_NO(T_g_>=1200,1)=Pn5(35.9916912,0.957170,-0.148032, 0.009974,-3.004088,T_g_(T_g_>=1200))/M_NO;
    
    c_N2O(T_g_<1400,1)=Pn5(27.67,51.148,-30.64, 6.847911, -0.157906,T_g_(T_g_<1400))/M_N2O;
    c_N2O(T_g_>=1400,1)=Pn5(60.30274,1.034566,-0.192997, 0.012540, -6.860254,T_g_(T_g_>=1400))/M_N2O;
    
    c_H2O(T_g_<1700,1)=Pn5(30.092,6.832514,6.793435,-2.534480,0.082139,T_g_(T_g_<1700))/M_H2O;
    c_H2O(T_g_>=1700,1)=Pn5(41.96426,8.622053,-1.499780,0.098119,-11.15764,T_g_(T_g_>=1700))/M_H2O;
    c_CH2O(T_g_<1200,1)=Pn5(5.19,93.2,-44.85, 7.882279,0.551175,T_g_(T_g_<1200))/M_CH2O;
    c_CH2O(T_g_>=1200,1)=Pn5(71.35268,6.174497,-1.191090, 0.079564,-15.58917,T_g_(T_g_>=1200))/M_CH2O;
    c_CO(T_g_<1300,1)=Pn5(25.56759,6.096130,4.054656, -2.671301, 0.131021,T_g_(T_g_<1300))/M_CO;
    c_CO(T_g_>=1300,1)=Pn5(35.15070,1.300095,-0.205921, 0.013550, -3.282780,T_g_(T_g_>=1300))/M_CO;
    
    c_gHMX=Pn2(0.669,77.82,T_g_)/M_HMX;
    
    C_BF3(T_g_<1000,1) = Pn5(21.28631, 130.3006, -109.9919, 34.28838, -0.073386,T_g_(T_g_<1000))/ M_BF3;
    C_BF3(T_g_>=1000,1) =  Pn5(81.23696, 1.096330, -0.226830, 0.015981, -6.366625,T_g_(T_g_>=1000))/M_BF3;
    
    c_gTf(T_g_<1100,1)=Pn5(43.55126,175.9079,-138.7331,40.35619,-0.381260,T_g_(T_g_<1100))/M_Tf;
    c_gTf(T_g_>=1100,1)=Pn5(129.9776,1.690918,-0.340087,0.023448,-10.83204,T_g_(T_g_>=1100))/M_Tf;
    
    C_CF2(T_g_<600,1) = Pn5(8.825772, 125.3652, -129.4940, 50.33101, 0.259749,T_g_(T_g_<600))/ M_CF2;
    C_CF2(T_g_>=600,1) =  Pn5(59.34753, -2.317176, 0.890518, -0.055879, -3.467545,T_g_(T_g_>=600))/M_CF2;
    
    C_C= Pn5(21.17510, -0.812428, 0.448537, -0.043256, -0.013103,T_g_)/M_C;
    
    C_C2(T_g_<700,1) =  Pn5(123.9092,-348.0067,485.0971, -232.7994, -1.240298,T_g_(T_g_<700))/M_C2;
    C_C2(T_g_>=700,1) =  Pn5(30.50408,5.445811,-0.853373, 0.065641, 0.814750,T_g_(T_g_>=700))/M_C2;
    
    C_HMXproducts=c_NO2*Y_NO2_f+c_NO*Y_NO_f+c_N2O*Y_N2O_f+c_H2O*Y_H2O_f+c_CO*Y_CO_f+c_CH2O*Y_CH2O_f;
    cg=C_HMXproducts.*Y_HMXprod_+c_gHMX.*Y_gHMX_+c_gTf.*Y_gTf+C_BF3.*Y_BF3 + C_CF2.*Y_CF2 + C_C.*Y_C+C_C2.*Y_C2;
    % теплоемкости к-веществ
    
%     c_cTf(T_g_<TmTf,1)=1.04;
%     c_cTf(T_g_>=TmTf,1)= (0.61488+0.001194*T_g_);
   c_cTf= (0.61488+0.001194*T_g_);
    c_B = (10.18574 + 29.24415*tau(T_g_) -18.02137*tau(T_g_).^2 +4.212326*tau(T_g_).^3 )/10.8 ;
    cc=c_cTf.*Y_cTf+ c_B.*Y_B;
    c_gf=((1-fg).*roc.*cc+fg.*rog.*cg)./ro_gf;
    
    % теплопроводности
    lc=2.5e-3*Y_cTf  + 0.27*Y_B;
    l_HMXproducts = 2e-4*Y_NO2_f+2.59e-4*Y_NO_f+1.74e-4*Y_N2O_f+2.79e-4*Y_H2O_f+2.5e-4*Y_CO_f+ 2e-4*Y_CH2O_f; % l_BF3
    lg=l_HMXproducts*Y_HMXprod_+12.5e-4*Y_gHMX_+2.5e-4*Y_gTf +2e-4*(Y_BF3  +Y_C+Y_C2 +Y_CF2);
    l_gf=((1-fg).*lc.*cc+fg.*lg.*cg)./((1-fg).*cc+fg.*cg);
    
    % кинетмка
    %Tf->C2F4
    % C2F4->2CF2
    % B+3CF2->2BF3+3C
    % 2C->C2
     selectHMX=1;
    switch selectHMX
        case 1
            G_gHMX =fg.*ro_gf.*Y_gHMX_.*Ar(10^14.2, 165268)/M_HMX; %
        case 2
            G_gHMX =fg.*ro_gf.*Y_gHMX_.*Ar(10^12.5, 158992)/M_HMX;
        case 3
            G_gHMX =fg.*ro_gf.*Y_gHMX_.*Ar(5*10^13, 158992)/M_HMX;
            
    end
    G_Tf=(1-fg).*ro_cTf.*Y_cTf.*Ar(kTf,ETf)/M_Tf;
    G_CF2dec =fg.*ro_gf.*Y_gTf.*T_g_.^(-0.5).*Ar(7.82*10^15,2.33*10^5)/M_C2F4;
    G_B_CF2 = fg*ro_gf.*Y_CF2.*(Y_B*ro_B).^3.*T_g_.^(0.5).*Ar(4*10^14,8.37*10^4)/(M_CF2*M_B^3);
    G_C =fg.*ro_gf.^2.*Y_C.^2.*T_g_.^(-1.6).*Ar(1.80*10^21,0)/M_C;
    
    w_gHMX=  -M_HMX*G_gHMX;
    w_gHMXprod = M_HMX*G_gHMX;
    %Tf->C2F4
    w_cTf=-M_Tf*G_Tf;
    w_gTf = -M_Tf*w_cTf;
    % B+3CF2->BF3+3C
    w_B=-M_B*G_B_CF2;
    % C2F4->2CF2, B+3CF2->2BF3+3C
    w_CF2=2*M_CF2*G_CF2dec-3*M_CF2*G_B_CF2;
    % B+3CF2->2BF3+3C
    w_BF3 =2*M_BF3* G_B_CF2;
    % B+3CF2->2BF3+3C, 2C->C2
    w_C = -2*M_C*G_C +3*M_C*G_B_CF2;
    w_C2 = M_C2*G_C;
    
    G_condphase= -w_cTf-w_B;
    % тепловыделение
    Q_react_g= Q_gHMX*w_gHMX/(wm*dHm)+ ( -w_cTf*824/M_Tf - w_gTf*658.56/M_Tf+ ...       $ -258
        + w_CF2*dH_CF2/M_CF2 + w_BF3*dH_BF3/M_BF3  + dH_C*w_C/M_C +dH_C2*w_C2/M_C2)*1e3/(wm*dHm);
    %disp(Q_react_g);
    Pe_g=cfm*m.*l./l_gf;
    Pe_dk = m*l./(D*rog);
    res_T_g_ = zeros(n,1);
    res_Y_gHMX_ = zeros(n,1);
    res_Y_HMXprod_=zeros(n,1);
 % EQUATIONS_____________________________   
    ind = 2:n-1;
    res_T_g_(1) = T_g_(1)-T_fout;
    res_T_g_(ind) =  (Pe_gm./Pe_g(ind)).*(T_g_(ind+1)-2*T_g_(ind)+T_g_(ind-1))/dx^2 + 0*(1./Pe_dk(ind)).*(c_gf(ind)/cfm).*ro_gf(ind).*((T_g_(ind+1)-T_g_(ind-1))/(2*dx)).*Y_gHMX_(ind).*(Y_gHMX_(ind+1)-Y_gHMX_(ind-1))/(2*dx) - ...
        (Pe_gm)^0.5 .* (c_gf(ind)/cfm).*(T_g_(ind+1)-T_g_(ind))/(2*dx) -  Q_react_g(ind);   

    res_T_g_(n) = ( T_g_(n)-T_g_(n-2))/(2*dx) ;
    
    res_Y_gHMX_(1) = Y_gHMX_(1)-Y_gHMXout;    
    res_Y_gHMX_(ind) =(Pe_gm./Pe_dk(ind)).*(Y_gHMX_(ind+1)-2*Y_gHMX_(ind)+Y_gHMX_(ind-1))/dx^2- (Pe_gm)^0.5 .*(Y_gHMX_(ind+1)-Y_gHMX_(ind-1))/(2*dx)+w_gHMX(ind)/wm; %w_gHMX  Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
    res_Y_gHMX_(n) =( Y_gHMX_(n)-Y_gHMX_(n-2))/(2*dx);
    
    res_Y_HMXprod_(1) = Y_HMXprod_(1)-Y_HMXprod0;
    res_Y_HMXprod_(ind) =(Pe_gm./Pe_dk(ind)).*(Y_HMXprod_(ind+1)-2*Y_HMXprod_(ind)+Y_HMXprod_(ind-1))/dx^2- (Pe_gm)^0.5 .*(Y_HMXprod_(ind+1)-Y_HMXprod_(ind-1))/(2*dx)+w_gHMXprod(ind)/wm; %w_gHMX  Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
    res_Y_HMXprod_(n) =( Y_HMXprod_(n)-Y_HMXprod_(n-2))/(2*dx);
    
    res = [res_T_g_;res_Y_gHMX_; res_Y_HMXprod_;];
    end
% SOLVERS__________________________________________________________________
% SOLVER 1
function [sol, error] = newtonraphson(F, x0, maxiter)
% Newton-Raphson solver for a system of nonlinear equations, where F is the
% handle to the system of equations and x0 is the starting point for the
% solver. The output of F is a column vector, x0 is a column vector.
% Maxiter species the maximum number of iterations, default is 1000.

% initialise values
error1 = 1e8;
x = x0;
iter = 0;
if nargin < 3
    maxiter = 1e3;
end
error = zeros(maxiter, 1);

sol = F(x0); % compute test solution
nvar = length(sol); % compute size of the system
h = 1e-4 .* ones(nvar, 1); % size of step for derivative computation

while error1 > 1e-12
    iter = iter+1; % update iteration
    f = F(x); % evaluate function at current point
    
    J = jacobiannum(F, x, h); % compute the Jacobian
    y = -J\f;  % solve the linear equations
    x = x + y; % move the solution
    
    % calculate errors
    error1 = sqrt(sum(y.^2));
    error(iter) = sqrt(sum(f.^2));
    
    % break computation if maximum number of iterations is exceeded
    if iter == maxiter
        warning('Maximum number of iterations (%i) exceeded, solver may not have converged.', maxiter);
        break;
    end 
end

% return solution and error for each iteration
sol = x;
error = error(1:iter);

end

function J = jacobiannum(F,x,h)
% Computes the Jacobian matrix of the function F at the point x, where h is
% the step size to take on each dimension (has to be small enough for
% precision). Note that the result of f and vectors x and h must be column
% vectors of the same dimensions.

J = (F(repmat(x,size(x'))+diag(h))-F(repmat(x,size(x'))))./h';

end
% SOLVER 2
function [f,xhist] = newton(fun,x0,varargin)
% Solves a system of nonlinear equations  fun(x) = 0 by Newton's method 
% (a.k.a. the Newton-Raphson method) Only real-valued solutions accepted.
% Usage:
%  x = newton(fun,x0)
%    fun is a handle to a function returning a column vector f of function 
%    values and optionally the Jacobian matrix J, where J(i,j) = df(i)/dx(j)
%       f = fun(x) or [f,J] = fun(x) 
%    x0 is the column vector of starting values for the search.
%    x is the solution, if one is found. Otherwise, newton issues a 
%       warning; "No convergence" and returns a vector of NaNs
%    f, x0 and x are all column vectors of length n. J is n by n.
%  x = newton(fun,x0,name,value, ...)
%    allows the user to control the iteration.  
%    The following names are allowed:
%      bounds: n by 2 matrix where bounds(i,1) and bounds(i,2) are 
%        lower and upper bounds, respectively, for variable i.   
%        Use -Inf and/or Inf to indicate no bound.  Default: No bounds.
%      maxiter: Maximum number of iterations,  Default: 50
%      tolx:    Minimum length of last iteration step. Default: 1e-4
%      tolfun:  Minimum value for norm(f). Default: 1e-5
%    Example: x = newton(fun,x0,'maxiter',100,'tolfun',1e-3) 
%  [x,xhist] = newton(fun,x0)
%     xhist contains the search history. x(k,i) is x(i) at iteration k
%   See also: - broyden 
%               (https://se.mathworks.com/matlabcentral/fileexchange/54667)
%             - fsolve (optimization toolbox)
% Are Mjaavatten, November 2020
% Version 2.0: January 2021
%   Option for numeric Jacobian
%   Option for bounds and for modifying iteration parameters
  
  % Parse any name-value pairs that will override the defaults:
  p = inputParser;
  addParameter(p,'tolfun',1e-5);  % Last argument is the default value
  addParameter(p,'tolx',1e-4);
  addParameter(p,'maxiter',50);
  addParameter(p,'bounds',[]);
  parse(p,varargin{:})
  opt = p.Results;    
  % Fields of opt take the default values unless overriden by a 
  % name/value pair
  
  bnds = ~isempty(opt.bounds);  % Bounds on x
  if any(diff(opt.bounds')<0)
    error('newton:bounds',...
      'bounds(i,2)-bounds(i,1) must be > 0 for all i')
  end  
  if bnds && ~all(isreal(opt.bounds))
    error('newton:complexbounds','Bounds cannot be complex numbers')
  end
  no_jacobian = nargout(fun) < 2;
  x0 = x0(:);  % Make sure x0 is a column vector
  if no_jacobian
    f = fun(x0);
    J = jacobi(fun,x0);
  else
    [f,J] = fun(x0);
  end
  x = x0(:);  % Make sure x0 is a column vector
  if ~(size(x) == size(f))
    error('newton:dimension',...
      'fun must return a column vector of the same size as x0')
  end  
  xhist = zeros(opt.maxiter,length(x));
  for i = 1:opt.maxiter
    xhist(i,:) = x';  % Search history
    if any(isnan(J(:))) || rcond(J) < 1e-15
      warning('newton:singular',...
        'Singular jacobian at iteration %d.\n Iteration stopped.',i);
      x = NaN*x0;
      return
    end     
    dx = -J\f(:);
    if bnds
      xnew = real(x + dx);
      xnew = min(max(xnew,opt.bounds(:,1)),opt.bounds(:,2));
      dx = xnew-x;
    end
    x = x + dx; 
    if norm(f(:))<opt.tolfun && norm(dx)<opt.tolx
      xhist = xhist(1:i,:);
      if norm(fun(real(x)))> opt.tolfun  % Complex part not negligible
        warning('newton:complex','Converged to complex solution.')
      else
        x = real(x);
      end
      return
    end    
    if no_jacobian
      f = fun(x);
      J = jacobi(fun,x);
    else
      [f,J] = fun(x);
    end
  end
  if any(abs(xhist(end-1,:)-xhist(end,:))<= 1e-8) 
    s1 = 'Search stuck. ';
    s2 = 'Possibly because Newton step points out of bounded region.';
    if bnds
      warning('newton:stuck',[s1,s2,'\nNo convergence']);
    else
      warning('newton:stuck',[s1,'\nNo convergence']);
    end
    fprintf('Last x:\n')
    for k = 1:numel(x)
      fprintf('%f\n',x(k))
    end
    x = x*NaN;
  else
    warning('newton:noconv','No convergence.')
    x = x*NaN;
  end
end
function J = jacobi(f,x,delta)
% Simple Jacobian matrix estimate using central differences
%  f: function handle
%  x: column vector of independet variables
%  delta: Optional step leghth.  Default: 1e-7*sqrt(norm(x))
%
% See also: John D'Errico's jacobianest.m in
%   https://www.mathworks.com/matlabcentral/fileexchange/13490
%   (Robust and high accuracy)
  if nargin < 3
      delta = 1e-7*sqrt(norm(x));
  end
  y0 = feval(f,x);
  n = length(y0);
  m = length(x);
  J = zeros(n,m);
  X = repmat(x,[1,m]);
  d = delta/2*eye(m);
  for i = 1:m
    J(:,i) = (f(X(:,i)+d(:,i))-f(X(:,i)-d(:,i)))/delta;
  end
end

end