problem in the for - loop

Here is the problem

Ive got a FEM code

function main1
clc
  n = 100;
  x0 = ones(3*n,1);
  sol = fsolve(@(x)fun(x,n),x0);
  norm(fun(sol,n))
  x = ((1:n)-1)/(n-1);
  plot(x,sol(1:n))
end
function res = fun(z,n)
  eta=1.0;
  beta = 1.0;
  x = ((1:n)-1)/(n-1);
  dx = 1/(n-1);
  
  y1 = z(1:n);
  y2 = z(n+1:2*n);
  y3 = z(2*n+1:3*n);
 for i=1:length(y1)
  Y1=y1(i);
 end
  for i=1:length(y2)
  Y2=y2(i);
 end
  res_y1 = zeros(n,1);
  res_y2 = zeros(n,1);
  res_y3 = zeros(n,1);
 
  res_y1(1) = y1(1)-10.0;
  for i=2:n-1
      
    res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) + exp(-5/Y2);
  end
  res_y1(n) = y1(n);
  res_y2(1) = y2(1);
  for i = 2:n-1
    res_y2(i) = (y2(i+1)-2*y2(i)+y2(i-1))/dx^2 - y1(i).^2;
  end
  res_y2(n) = y2(n)-eta;
  res_y3(1) = y3(1);
  for i=2:n-1
    res_y3(i) = (y3(i+1)-2*y3(i)+y3(i-1))/dx^2 - Y2;
  end
  res_y3(n) = y3(n)-1.0;
  res = [res_y1;res_y2;res_y3];
end

It seems that using

res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) + exp(-5/Y2);

and

res_y1(i) = (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) + exp(-5/y2(i));

must be eqiuvalent, but NO. The results differ

Answer 1

[Gleb]

Ok, thank you. I rewrote the code. It is more readable now. But "Unable to perform assignment because the left and right sides have a different number of elements.

Error in Gas_system_4v5/fun (line 219)

res_T_g_(ind) = (Pe_gm./Pe_g(ind)).*(T_g_(ind+1)-2*T_g_(ind)+T_g_(ind-1))/dx^2 - ...

"

again. The answer i assume is simple: cc, cg, Pe_g are matrixes n x n, not vectors, its because elementwise to n vectors gives matrix.

function Gas_system_4v5
clc
close all
n =50;
global T_fout Y_gHMXout Pe_g Pe_dk w_gHMX Pe_gm;
  x0 = [800*ones(n,1);1*ones(n,1);1*ones(n,1)];
  x = ((1:n)-1)/(n-1);
  Pe_gm=4.8157e+004;
 AlgoSelector=1;
 switch AlgoSelector
     case 1
 Algochise='levenberg-marquardt'; %  'trust-region-dogleg'
     case 2
 Algochise='trust-region-dogleg' ;
     case 3
 Algochise='trust-region' ;
 end
%  options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-5, 'FiniteDifferenceStepSize', 10^-6);
 % options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-3, 'FiniteDifferenceStepSize', 10^-4);
   % 'OptimalityTolerance', 10^-11,
  %options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'FunctionTolerance', 10^-4, 'OptimalityTolerance', 10^-3, 'StepTolerance', 10^-5  );    % optimoptions
  options=optimset( 'Display','iter', 'MaxFunEvals', 10^8, 'MaxIter', 10^6, 'Algorithm', Algochise, 'TolFun', 10^-7);   
  [sol,fval] = fsolve(@(x)fun(x,n),x0,  options);
  disp(fval);
  norm(fun(sol,n))
   figure
  plot(x*(  Pe_gm)^-0.5,sol(1:n))
  hold on
  figure
    plot(x,sol(n+1:2*n))
    hold on
     figure
    plot(x,sol(2*n+1:3*n))
    hold on
disp(Pe_g);
disp(Pe_dk);
function res = fun(Yg,n)
  x = ((1:n)-1)/(n-1);
  dx = 1/(n-1);
   T_fout =750;
Y_gHMXout =0.8;
 cfm=5*10^4;
  T_g_ = Yg(1:n);
  Y_gHMX_ = Yg(n+1:2*n);
  Y_HMXprod_ = Yg(2*n+1:3*n);
%   for i=1:length(T_g_)
% T_g= T_g_(i);
%   end
%    %for i=1:n
%   for i=1:length(Y_gHMX_) 
% Y_gHMX= Y_gHMX_(i);
%   end
%global T_fout  L_g Y_gHMXout  Y_B0 x  Q_react_g
% format long e
Y_HMXprod0=0.2;
f1Y_cTfout=0.0;
Y_B0 = 0.0;
 f_gout=0.99;
l=0.5;
 f1Y_B0=Y_B0 * (1-f_gout);
 selectHMX=1;
 Choice1=1;
 solverselect=1;
 Selector=2 ;
 
MAX = 100;
bcinterval=[-10^2 10^2];   
     
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Y_HMXprod=0.3;
Y_gTf=0;
Y_BF3 =0;
Y_CF2=0;
Y_C=0;
Y_C2=0;
Y_cTf=0.01;
 Y_B=0.01;
 fg=0.7;
 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;
T_g=T_g_;
Y_HMXprod=Y_HMXprod_;
Y_gHMX=Y_gHMX_;
  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)=Pn5(16.108,75.89,-54.3874,14.30777,0.239423, T_g_(T_g_<1200))/M_NO2;
c_NO2(T_g_>=1200)=Pn5(56.82541,0.738053,-0.144721, 0.009777,-5.459911,T_g_(T_g_>=1200))/M_NO2;
c_NO(T_g_<1200)=Pn5(23.83, 12.589, -1.139, -1.497459, 0.214194, T_g_(T_g_<1200))/M_NO ;
c_NO(T_g_>=1200)=Pn5(35.9916912,0.957170,-0.148032, 0.009974,-3.004088, T_g_(T_g_>=1200))/M_NO  ;  
c_N2O(T_g_<1400)=(Pn5(27.67,51.148,-30.64, 6.847911, -0.157906, T_g_(T_g_<1400))/M_N2O );
c_N2O(T_g_>=1400)=Pn5(60.30274,1.034566,-0.192997, 0.012540, -6.860254,  T_g_(T_g_>=1400))/M_N2O  ;
c_H2O(T_g_<1700)=Pn5(30.092,6.832514,6.793435,-2.534480,0.082139,T_g_(T_g_<1700) )/M_H2O;
c_H2O(T_g_>=1700)=Pn5(41.96426,8.622053,-1.499780,0.098119,-11.15764, T_g_(T_g_>=1700) )/M_H2O;
c_CH2O(T_g_<1200)=Pn5(5.19,93.2,-44.85, 7.882279,0.551175, T_g_(T_g_<1200))/M_CH2O ;
c_CH2O(T_g_>=1200)=Pn5(71.35268,6.174497,-1.191090, 0.079564,-15.58917,T_g_(T_g_>=1200) )/M_CH2O ;
c_CO(T_g_<1300)=Pn5(25.56759,6.096130,4.054656, -2.671301, 0.131021, T_g_(T_g_<1300))/M_CO ;
c_CO(T_g_>=1300)=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) = Pn5(21.28631, 130.3006, -109.9919, 34.28838, -0.073386, T_g_(T_g_<1000))/ M_BF3;
C_BF3(T_g_>=1000) =  Pn5(81.23696, 1.096330, -0.226830, 0.015981, -6.366625, T_g_(T_g_>=1000))/M_BF3;
c_gTf(T_g_<1100)=Pn5(43.55126,175.9079,-138.7331,40.35619,-0.381260,  T_g_(T_g_<1100))/M_Tf;
c_gTf(T_g_>=1100)=Pn5(129.9776,1.690918,-0.340087,0.023448,-10.83204,  T_g_(T_g_>=1100))/M_Tf;
C_CF2(T_g_<600) =Pn5(8.825772, 125.3652, -129.4940, 50.33101, 0.259749, T_g_(T_g_<600))/ M_CF2;
C_CF2(T_g_>=600) =  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) = Pn5(123.9092,-348.0067,485.0971, -232.7994, -1.240298, T_g_(T_g_<700))/M_C2;
C_C2(T_g_>=700) = 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.04;
c_cTf(T_g_>=TmTf)= (0.61488+0.001194*1000*tau(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
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); %0*( -w_cTf*824/M_Tf - w_gTf*658.56/M_Tf + 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);
  %________________________________________________________________________
  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 - ...
    (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-1))/dx ;
  %________________________________________________________________________
    ind = 2:n-1;
  res_Y_gHMX_(1) = Y_gHMX_(1)-Y_gHMXout;

    res_Y_gHMX_(ind) =(1./Pe_dk(ind))*(Y_gHMX_(ind+1)-2*Y_gHMX_(ind)+Y_gHMX_(ind-1))/dx^2-(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-1))/dx;
  %________________________________________________________________________
  ind = 2:n-1;
  res_Y_HMXprod_(1) = Y_HMXprod_(1)-Y_HMXprod0;     

 res_Y_HMXprod_(ind) =(1./Pe_dk(ind))*(Y_HMXprod_(ind+1)-2*Y_HMXprod_(ind)+Y_HMXprod_(ind-1))/dx^2-(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-1))/dx;
  
  res = [res_T_g_;res_Y_gHMX_; res_Y_HMXprod_;];
end
end

[Cris LaPierre]

You are correct. i had to make sure to use the elementwise operator (.* and ./) in calculating those variables. I also had to be sure all vectors were nx1 to avoid this issue, so I assign using c_NO2(T_g_<1200,1)=...

Here's what my solution looks like. Notice I also avoid using global variables. I unfortunately have no idea if the result is correct.

Gas_system_4v3

                                        First-Order                    Norm of 
 Iteration  Func-count    Residual       optimality      Lambda           step
     0         151     9.18573e+06        7.38e+05         0.01
     1         302     8.51224e+06        1.17e+07        0.001         63.702
     2         453          273706        1.38e+06       0.0001         62.519
     3         604         12191.6        8.26e+04        1e-05        56.9203
     4         755         303.796        1.67e+03        1e-06        53.7196
     5         906         20.3801             475        1e-07        201.989
     6        1057       0.0823557              27        1e-08        53.7634
     7        1208     5.60334e-06          0.0221        1e-09        5.96644
     8        1359     2.08205e-13        1.53e-06        1e-10     0.00192897

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

ans = 4.5629e-07

function Gas_system_4v3
n =50;

x0 = [800*ones(n,1);1*ones(n,1);1*ones(n,1)];
x = ((1:n)-1)/(n-1);
Pe_gm=4.8157e+004;
AlgoSelector=1;
switch AlgoSelector
    case 1
        Algochise='levenberg-marquardt'; %  'trust-region-dogleg'
    case 2
        Algochise='trust-region-dogleg' ;
    case 3
        Algochise='trust-region' ;
end
%  options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-5, 'FiniteDifferenceStepSize', 10^-6);
% options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-3, 'FiniteDifferenceStepSize', 10^-4);
% 'OptimalityTolerance', 10^-11,
%options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'FunctionTolerance', 10^-4, 'OptimalityTolerance', 10^-3, 'StepTolerance', 10^-5  );    % optimoptions
options=optimset( 'Display','iter', 'MaxFunEvals', 10^8, 'MaxIter', 10^6, 'Algorithm', Algochise, 'TolFun', 10^-7);
[sol,fval] = fsolve(@(x)fun(x,n),x0,  options);
% Call the function with the solution values to capture Pe_g and Pe_dk
[~,Pe_g, Pe_dk] = fun(sol,n);
% disp(fval);
norm(fun(sol,n))

figure
plot(x*(Pe_gm)^-0.5,sol(1:n))
hold on

figure
plot(x,sol(n+1:2*n))
hold on

figure
plot(x,sol(2*n+1:3*n))
hold on

% disp(Pe_g);
% disp(Pe_dk);

    function [res, Pe_g, Pe_dk] = fun(Yg,n)
    x = ((1:n)-1)/(n-1);
    dx = 1/(n-1);
    T_fout =750;
    Y_gHMXout =0.8;
    cfm=5*10^4;
    T_g_ = Yg(1:n);
    Y_gHMX_ = Yg(n+1:2*n);
    Y_HMXprod_ = Yg(2*n+1:3*n);
    %   for i=1:length(T_g_)
    % T_g= T_g_(i);
    %   end
    %    %for i=1:n
    %   for i=1:length(Y_gHMX_)
    % Y_gHMX= Y_gHMX_(i);
    %   end

    % format long e
    
    Y_HMXprod0=0.2;
    f1Y_cTfout=0.0;
    Y_B0 = 0.0;
    f_gout=0.99;
    l=0.5;
    
    f1Y_B0=Y_B0 * (1-f_gout);
    
    selectHMX=1;
    Choice1=1;
    solverselect=1;
    Selector=2 ;
    
    MAX = 100;
    bcinterval=[-10^2 10^2];
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Y_HMXprod=0.3;
    Y_gTf=0;
    Y_BF3 =0;
    Y_CF2=0;
    Y_C=0;
    Y_C2=0;
    Y_cTf=0;
    Y_B=0;
    fg=0.7;
    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;
    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_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
    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)+0*( -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);
    
    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);   %  Q_gHMX*Y_gHMX_(i).*ro_gf.*Ar(10^14.2, 165268) 1*

    res_T_g_(n) = ( T_g_(n)-T_g_(n-1))/dx ;
    res_Y_gHMX_(1) = Y_gHMX_(1)-Y_gHMXout;    
    res_Y_gHMX_(ind) =(1./Pe_dk(ind)).*(Y_gHMX_(ind+1)-2*Y_gHMX_(ind)+Y_gHMX_(ind-1))/dx^2-(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-1))/dx;
    
    
    res_Y_HMXprod_(1) = Y_HMXprod_(1)-Y_HMXprod0;
    %  w_gHMXprod=  fg*Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ));
    res_Y_HMXprod_(ind) =(1./Pe_dk(ind)).*(Y_HMXprod_(ind+1)-2*Y_HMXprod_(ind)+Y_HMXprod_(ind-1))/dx^2-(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-1))/dx;
    
    res = [res_T_g_;res_Y_gHMX_; res_Y_HMXprod_;];
    end
end

[Gleb]

Yap, it works!

But i have a question. It seems that addying 1 to the structure like C_CF2(T_g_<600,1) you transpose the vector. Isn,t it?

and what the meaning of [~,Pe_g, Pe_dk] = fun(sol,n), why this functions are calling at the begining of code?

[Cris LaPierre]

I put this in the answer: "I also had to be sure all vectors were nx1 to avoid this issue, so I assign using c_NO2(T_g_<1200,1)=..."

This syntax places the returned values into a single column instead of a row. This was necessary to prevent implicit explansion turning the results of some calculations into 50x50 matrices.

As for the second question, this was how I elected to make the 2 variables returned by the function (Pe_g, Pe_dk) available for the disp commands (I commented them out). You used a global variable, but it's best to avoid globals if at all possible.

If you remove the disp commands, you can remove [~,Pe_g, Pe_dk] = fun(sol,n) as well.

[Gleb]

Thank you. Problem solved. Regards.

Answer 2

[Cris LaPierre]

Why would you expect them to be the same?

When you get around to calculation res_y1(i), Y2 has the value of y2(length(y1)). If you substitute that in, you'll see the last part fo the equation is very different

exp(-5/y2(length(y1)));
%vs
exp(-5/y2(i));

The first uses the same value of Y2 for every loop, while the second uses a different element of y2 for each loop.

[Gleb]

I corrected the code. I write

  for i=1:length(y2)
  Y2=y2(i);
  end
 

but the problem remains. I just want some functions of unknown variables to be calculated and then substituted into equations because these functions can be too huge to fit directly into the equations themselves. for example

    (y1(i+1)-2*y1(i)+y1(i-1))/dx^2 + (y1(i+1)-y1(i-1))/(2*dx) +Q

where

    Q =f(y(i), y(i), y(3))

[Cris LaPierre]

The underlying problem appears to still be the same. In one case, Y2 is always the same value for every loop while in the other, y2(i) will be a different value in every loop.

[Gleb]

so the only way is to write Q inside the for-loop? but Q in my problems often use some

function

and functions can not be used inside the loop

[Cris LaPierre]

You can absolutely use a function inside a for loop. You can't define an in-file function in the loop, but that's a different issue.

Depending what Q is, you could use a in-file function, or perhaps an anonymous function.

[Gleb]

Ok, thank you, i placed the anon. function inside the for-loop of the the 1st equation and it works, but i discovered that functions inside the 1st loop are not available for the 2nd for-loop.

Example

  res_Y1(1) = Y1(1)-Y0;
  for i = 2:n-1
      Q=@(...)(...) % some function of Y1(i) Y2(i) and parameters
 res_Y1(i) =(Y1(i+1)-2*Y1(i)+Y1(i-1))/dx^2-Q
  end
  res_Y1(n) =( Y1(n)-Y1(n-1))/dx;
 
  res_Y2(1) = Y2(1)-Y20;
 for i = 2:n-1
    
 res_Y2(i) =(Y2(i+1)-2*Y2(i)+Y2(i-1))/dx^2+Q
  end
 res_Y2(n) =(Y2(n)-Y2(n-1))/dx;
  
  res = [ res_Y1 ;  res_Y2];

so it seems to me that the loop for the 2nd equation "thinks" Q is a constant

[Cris LaPierre]

I guess it's all in how you define your anonymous function. If you define your anonymous function in a variable, you can use it as a function elsewhere in your code.

% define an anonymous function
sqr = @(x) x.^2;
% use the anonymous function
a = sqr(5)

a = 25

[Gleb]

So is it correct in my example to use such structure

  res_Y1(1) = Y1(1)-Y0;
  for i = 2:n-1
A=@(Y1(i), Y2(i))(...) % some expressions
B=@(Y1(i), Y2(i))(...)
Q=@(Y1(i), Y2(i))(A(Y1(i), Y2(i))+B(Y1(i), Y2(i))) 
 res_Y1(i) =(Y1(i+1)-2*Y1(i)+Y1(i-1))/dx^2-Q(Y1(i), Y2(i))
   end
  res_Y1(n) =( Y1(n)-Y1(n-1))/dx;
 % loop 2
   res_Y2(1) = Y2(1)-Y20;
 for i = 2:n-1
     res_Y2(i) =(Y2(i+1)-2*Y2(i)+Y2(i-1))/dx^2+Q(Y1(i), Y2(i))
 end
 ...

and have a hope that Q will not be a constant in other loops

[Cris LaPierre]

No, I don't think I would code it like this. Of course, you have left out a lot of the details, so it's hard to say for certain. Generally, I wouldn't define the anonymous function inside a for loop.

Also, you are not following the syntax shared in the example. You should indicate a variable used in the equation right after the @, then use that variable in the function. When you call the function, that is when you use your actual indexed variable.

So based on your example, I think it would look like this (untested).

res_Y1(1) = Y1(1)-Y0;
A=@(myY1, myY2) (yY1+myY); % some expressions
B=@(myY1, myY2) (yY1.*myY);
Q=@(myY1, myY2) A(myY1, myY2)+B(myY1, myY2);
for i = 2:n-1
    res_Y1(i) =(Y1(i+1)-2*Y1(i)+Y1(i-1))/dx^2-Q(Y1(i), Y2(i))
end
res_Y1(n) =( Y1(n)-Y1(n-1))/dx;
% loop 2
res_Y2(1) = Y2(1)-Y20;
for i = 2:n-1
    res_Y2(i) =(Y2(i+1)-2*Y2(i)+Y2(i-1))/dx^2+Q(Y1(i), Y2(i))
end
...

[Gleb]

Ok, thank you.

And it will be i think the last question. I have some functions (heat capacity) set at intervals. example

if myY1<1000 
A=@(myY1, myY2) (myY1+myY2^2); % some expressions
else
A=@(myY1, myY2) (2*myY1+3*myY2^2)
end

obviously i can not use this code because it will be "Undefined function or variable". So polynomials should "change on signal" from the loop.

[Cris LaPierre]

The variable myY1 is a variable that only exists in the anonymous function. When you call A, that is when you pass in the actual value.

Also note that here you have only defined functions. Be sure to also use them to perform an actual calculation.

if Y1(i)<1000 
    % define your function
    A=@(myY1, myY2) (myY1+myY2^2); % some expressions
    % use your function
    res(i) = A(Y1(i),Y2(i));
else
    % define your function
    A=@(myY1, myY2) (2*myY1+3*myY2^2)
    % use your function
    res(i) = A(Y1(i),Y2(i));
end

If you are finding functions confusing, perhaps Function Basics will be helpful.

[Gleb]

But in your code Y1(i) is used before the loop is started.

[Cris LaPierre]

For that matter, Y1(i) is never defined in the code you shared. This is the problem with trying to simplify your example. If you want an answer that is guaranteed to work, you need to share your actual code.

[Gleb]

Here it is. since it is quite large, I used the "toy" examples in previous questions

On line 156 there is a polynom that should be "switched" on certain value of variable

function Gas_system_4v3
clc
close all
n =50;
global T_fout Y_gHMXout Pe_g Pe_dk w_gHMX;

T_fout =750;
Y_gHMXout =0.8;
  % x0 = ones(2*n,1);
 % x0 = [T_fout + 150*((1:n)-1)/(n-1);Y_gHMXout*((((1:n)-1)/(n-1)) -1)];
  x0 = [800*ones(n,1);1*ones(n,1);1*ones(n,1)];
  x = ((1:n)-1)/(n-1);
  Pe_gm=4.8157e+004;
 AlgoSelector=1;
 switch AlgoSelector
     case 1
 Algochise='levenberg-marquardt'; %  'trust-region-dogleg'
     case 2
 Algochise='trust-region-dogleg' ;
     case 3
 Algochise='trust-region' ;
 end
%  options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-5, 'FiniteDifferenceStepSize', 10^-6);
 % options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-3, 'FiniteDifferenceStepSize', 10^-4);
   % 'OptimalityTolerance', 10^-11,
  %options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'FunctionTolerance', 10^-4, 'OptimalityTolerance', 10^-3, 'StepTolerance', 10^-5  );    % optimoptions
  options=optimset( 'Display','iter', 'MaxFunEvals', 10^8, 'MaxIter', 10^6, 'Algorithm', Algochise, 'TolFun', 10^-7);   
  [sol,fval] = fsolve(@(x)fun(x,n),x0,  options);
  disp(fval);
  norm(fun(sol,n))
   figure
  plot(x*(  Pe_gm)^-0.5,sol(1:n))
  hold on
  figure
    plot(x,sol(n+1:2*n))
    hold on
     figure
    plot(x,sol(2*n+1:3*n))
    hold on

disp(Pe_g);
disp(Pe_dk);





function res = fun(Yg,n)
  x = ((1:n)-1)/(n-1);
  dx = 1/(n-1);
   T_fout =750;
Y_gHMXout =0.8;
 cfm=5*10^4;
  T_g_ = Yg(1:n);
  Y_gHMX_ = Yg(n+1:2*n);
  Y_HMXprod_ = Yg(2*n+1:3*n);
%   for i=1:length(T_g_)
% T_g= T_g_(i);
%   end
%    %for i=1:n
%   for i=1:length(Y_gHMX_) 
% Y_gHMX= Y_gHMX_(i);
%   end
%global T_fout  L_g Y_gHMXout  Y_B0 x  Q_react_g
% format long e

Y_HMXprod0=0.2;
f1Y_cTfout=0.0;
Y_B0 = 0.0;
 f_gout=0.99;
l=0.5;

 f1Y_B0=Y_B0 * (1-f_gout);

 selectHMX=1;
 Choice1=1;
 solverselect=1;
 Selector=2 ;
 
MAX = 100;
bcinterval=[-10^2 10^2];   
     
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Y_HMXprod=0.3;
Y_gTf=0;
Y_BF3 =0;
Y_CF2=0;
Y_C=0;
Y_C2=0;
Y_cTf=0;
 Y_B=0;
 fg=0.7;
 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;
  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+k2*tau.^2+k3*tau.^3+k4./tau(T_g).^2);
Ar=@(k,E, T_g)(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=@(T_g, Y_HMXprod, Y_gHMX)(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=@(T_g, Y_HMXprod, Y_gHMX)(Mg(T_g, Y_HMXprod, Y_gHMX).*P./(R.*T_g));
roc = ro_cTf*Y_cTf + ro_B*Y_B ;
ro_gf=@(T_g, Y_HMXprod, Y_gHMX)((1-fg)*roc+fg*rog(T_g, Y_HMXprod, Y_gHMX));

% теплоемкости газов

if T_g_(i)<1200
c_NO2=@(T_g)(Pn5(16.108,75.89,-54.3874,14.30777,0.239423)/M_NO2);

else
c_NO2=@(T_g)(Pn5(56.82541,0.738053,-0.144721, 0.009777,-5.459911)/M_NO2);

end
if T_g<1200
c_NO=@(T_g)(Pn5(23.83, 12.589, -1.139, -1.497459, 0.214194)/M_NO) ;
else
c_NO=@(T_g)(Pn5(35.9916912,0.957170,-0.148032, 0.009974,-3.004088)/M_NO)  ;  
end
if T_g<1400
c_N2O=@(T_g)(Pn5(27.67,51.148,-30.64, 6.847911, -0.157906)/M_N2O );
else
   c_N2O=@(T_g)(Pn5(60.30274,1.034566,-0.192997, 0.012540, -6.860254)/M_N2O ) ;
end
if T_g<1700
c_H2O=@(T_g)(Pn5(30.092,6.832514,6.793435,-2.534480,0.082139)/M_H2O);
else
    c_H2O=@(T_g)(Pn5(41.96426,8.622053,-1.499780,0.098119,-11.15764)/M_H2O);
end
if T_g<1200
c_CH2O=@(T_g)(Pn5(5.19,93.2,-44.85, 7.882279,0.551175)/M_CH2O) ;
else
    c_CH2O=@(T_g)(Pn5(71.35268,6.174497,-1.191090, 0.079564,-15.58917)/M_CH2O );
end
if T_g<1300
c_CO=@(T_g)(Pn5(25.56759,6.096130,4.054656, -2.671301, 0.131021)/M_CO) ;
else
c_CO=@(T_g)(Pn5(35.15070,1.300095,-0.205921, 0.013550, -3.282780)/M_CO) ;

end

c_gHMX=@(T_g)(Pn2(0.669,77.82)/M_HMX);
 if T_g<1000
C_BF3 = @(T_g)(Pn5(21.28631, 130.3006, -109.9919, 34.28838, -0.073386)/ M_BF3);
 else
C_BF3 =  @(T_g)(Pn5(81.23696, 1.096330, -0.226830, 0.015981, -6.366625)/M_BF3);
 end
 
 if T_g<1100
 c_gTf=@(T_g)(Pn5(43.55126,175.9079,-138.7331,40.35619,-0.381260)/M_Tf);
 else
 c_gTf=@(T_g)(Pn5(129.9776,1.690918,-0.340087,0.023448,-10.83204)/M_Tf);
 end
 
  if T_g<600
C_CF2 = @(T_g)(Pn5(8.825772, 125.3652, -129.4940, 50.33101, 0.259749)/ M_CF2);
 else
C_CF2 =  @(T_g)(Pn5(59.34753, -2.317176, 0.890518, -0.055879, -3.467545)/M_CF2);
  end
  
C_C= @(T_g)(Pn5(21.17510, -0.812428, 0.448537, -0.043256, -0.013103)/M_C); 

if T_g<700
C_C2 =  @(T_g)(Pn5(123.9092,-348.0067,485.0971, -232.7994, -1.240298)/M_C2);
else
C_C2 =  @(T_g)(Pn5(30.50408,5.445811,-0.853373, 0.065641, 0.814750)/M_C2);
end
  

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=@(T_g, Y_HMXprod, Y_gHMX)(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);
 % теплоемкости к-веществ 


 if T_g1<TmTf
 c_cTf=1.04;
 else
     c_cTf= @(T_g)(0.61488+0.001194*T_g);
 end
 c_B = @(T_g)(10.18574 + 29.24415*tau(T_g) -18.02137*tau(T_g).^2 +4.212326*tau(T_g).^3 )/10.8 ;    
cc(T_g)=c_cTf(T_g).*Y_cTf+ c_B(T_g).*Y_B; 
c_gf=@(T_g, Y_HMXprod, Y_gHMX)(((1-fg).*roc(T_g, Y_HMXprod, Y_gHMX).*cc(T_g)+fg.*rog(T_g, Y_HMXprod, Y_gHMX).*cg(T_g, Y_HMXprod, Y_gHMX))./ro_gf(T_g, Y_HMXprod, Y_gHMX));
T_ad=Q_gHMX*ro_gf*Y_gHMXout/C_HMXproducts;
% теплопроводности
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=@(T_g, Y_HMXprod, Y_gHMX)(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= @(T_g, Y_HMXprod, Y_gHMX)( ((1-fg).*lc.*cc(T_g)+fg.*lg(T_g, Y_HMXprod, Y_gHMX).*cg(T_g, Y_HMXprod, Y_gHMX))./((1-fg).*cc(T_g, Y_HMXprod, Y_gHMX)+fg.*cg(T_g, Y_HMXprod, Y_gHMX)));



% кинетмка
%Tf->C2F4
% C2F4->2CF2
% B+3CF2->2BF3+3C
% 2C->C2
switch selectHMX
    case 1
G_gHMX =@(T_g, Y_gHMX)(fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_gHMX.*Ar(10^14.2, 165268, T_g)/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=  @(T_g, Y_gHMX)(-M_HMX*G_gHMX(T_g, Y_gHMX));
w_gHMXprod = @(T_g, Y_gHMX)(M_HMX*G_gHMX(T_g, Y_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= @(T_g, Y_gHMX)((Q_gHMX*w_gHMX(T_g, Y_gHMX)/(wm*dHm)+0*( -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=@(T_g, Y_gHMX)(cfm*m.*l/l_gf(T_g, Y_gHMX)); 
Pe_dk = @(T_g, Y_gHMX)(m*l./(D*rog(T_g, Y_gHMX)));  
  res_T_g_ = zeros(n,1);
  res_Y_gHMX_ = zeros(n,1);
  res_Y_HMXprod_=zeros(n,1);

  res_T_g_(1) = T_g_(1)-T_fout;
  for i=2:n-1
  res_T_g_(i) =  (Pe_gm./Pe_g)*(T_g_(i+1)-2*T_g_(i)+T_g_(i-1))/dx^2 + 0*(1./Pe_dk).*(c_gf/cfm)*ro_gf.*((T_g_(i+1)-T_g_(i-1))/(2*dx)).*Y_gHMX_(i)*(Y_gHMX_(i+1)-Y_gHMX_(i-1))/(2*dx) - ...
  (Pe_gm)^0.5 * (c_gf/cfm)*(T_g_(i+1)-T_g_(i))/(2*dx) -  Q_react_g(T_g_(i), Y_gHMX_(i));   %  Q_gHMX*Y_gHMX_(i).*ro_gf.*Ar(10^14.2, 165268) 1*
  end
  res_T_g_(n) = ( T_g_(n)-T_g_(n-1))/dx ;
  res_Y_gHMX_(1) = Y_gHMX_(1)-Y_gHMXout;
  for i = 2:n-1

    res_Y_gHMX_(i) =(1./Pe_dk)*(Y_gHMX_(i+1)-2*Y_gHMX_(i)+Y_gHMX_(i-1))/dx^2-(Y_gHMX_(i+1)-Y_gHMX_(i-1))/(2*dx)+w_gHMX(T_g_(i), Y_gHMX_(i))/wm; %w_gHMX  Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
  end
  res_Y_gHMX_(n) =( Y_gHMX_(n)-Y_gHMX_(n-1))/dx;
  
 
  res_Y_HMXprod_(1) = Y_HMXprod_(1)-Y_HMXprod0;
 for i = 2:n-1
    %  w_gHMXprod=  fg*Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ));
 res_Y_HMXprod_(i) =(1./Pe_dk)*(Y_HMXprod_(i+1)-2*Y_HMXprod_(i)+Y_HMXprod_(i-1))/dx^2-(Y_HMXprod_(i+1)-Y_HMXprod_(i-1))/(2*dx)+w_gHMXprod(T_g_(i), Y_gHMX_(i))/wm; %w_gHMX  Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
  end
 res_Y_HMXprod_(n) =( Y_HMXprod_(n)-Y_HMXprod_(n-1))/dx;
  
  res = [res_T_g_;res_Y_gHMX_; res_Y_HMXprod_;];
end
end

[Cris LaPierre]

I see a few things.

1. First, your if statement uses T_g_(i), but i has not been defined.
2. You are defining a new function, c_NO2. However, you try to use it as a variable.
3. c_NO2 calls Pn5 with 5 inputs, but you defined it as having 6.

I suggest reviewing the link I shared about function basics. In particular, how to call a function. There are several anonymous functions created in the code that are not called properly.

• tau (no input provided in Pn5)
• Pn2 (has 3 inputs, but you call it with 2)
• Pn5 (has 6 inputs, but you call it with 5)
• Ar (has 3 inputs, but you call it with 2)
• ro_gf (has 3 inputs, but you often use it without any inputs

One more thing to figure out is what is T_g. You use it in your if statements conditional, but have never defined it. It is a variable used in some of your anonymous functions, but it only exists within the function's scope.

After fixing these issues, I think my approach to what you have in line 156 would use logical indexing instead. This allows me to compute a value for all values that meet the condition at once, rather than having to loop through each value.

c_NO2(T_g_<1200)=Pn5(16.108,75.89,-54.3874,14.30777,0.239423,T_g_(T_g_<1200))/M_NO2;
c_NO2(T_g_>=1200)=Pn5(56.82541,0.738053,-0.144721, 0.009777,-5.459911,T_g_(T_g_>=1200))/M_NO2;

Notice that c_NO2 is a variable that captures the result of the calculation performed by the function Pn5.

[Gleb]

I think I have corrected the code as per your recommendations. T_g is a variable for use in anonymous functions. But an error occurs "Index exceeds array bounds

Error in Gas_system_4v3>@(T_g)(c_cTf(T_g).*Y_cTf+c_B(T_g).*Y_B) (line 193)

cc=@(T_g)(c_cTf(T_g).*Y_cTf+ c_B(T_g).*Y_B);". i think there was some kind of confusion between functions and variables.

function Gas_system_4v3
clc
close all
n =50;
global T_fout Y_gHMXout Pe_g Pe_dk w_gHMX;
  x0 = [800*ones(n,1);1*ones(n,1);1*ones(n,1)];
  x = ((1:n)-1)/(n-1);
  Pe_gm=4.8157e+004;
 AlgoSelector=1;
 switch AlgoSelector
     case 1
 Algochise='levenberg-marquardt'; %  'trust-region-dogleg'
     case 2
 Algochise='trust-region-dogleg' ;
     case 3
 Algochise='trust-region' ;
 end
%  options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-5, 'FiniteDifferenceStepSize', 10^-6);
 % options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-3, 'FiniteDifferenceStepSize', 10^-4);
   % 'OptimalityTolerance', 10^-11,
  %options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'FunctionTolerance', 10^-4, 'OptimalityTolerance', 10^-3, 'StepTolerance', 10^-5  );    % optimoptions
  options=optimset( 'Display','iter', 'MaxFunEvals', 10^8, 'MaxIter', 10^6, 'Algorithm', Algochise, 'TolFun', 10^-7);   
  [sol,fval] = fsolve(@(x)fun(x,n),x0,  options);
  disp(fval);
  norm(fun(sol,n))
   figure
  plot(x*(  Pe_gm)^-0.5,sol(1:n))
  hold on
  figure
    plot(x,sol(n+1:2*n))
    hold on
     figure
    plot(x,sol(2*n+1:3*n))
    hold on

disp(Pe_g);
disp(Pe_dk);





function res = fun(Yg,n)
  x = ((1:n)-1)/(n-1);
  dx = 1/(n-1);
   T_fout =750;
Y_gHMXout =0.8;
 cfm=5*10^4;
  T_g_ = Yg(1:n);
  Y_gHMX_ = Yg(n+1:2*n);
  Y_HMXprod_ = Yg(2*n+1:3*n);
%   for i=1:length(T_g_)
% T_g= T_g_(i);
%   end
%    %for i=1:n
%   for i=1:length(Y_gHMX_) 
% Y_gHMX= Y_gHMX_(i);
%   end
%global T_fout  L_g Y_gHMXout  Y_B0 x  Q_react_g
% format long e

Y_HMXprod0=0.2;
f1Y_cTfout=0.0;
Y_B0 = 0.0;
 f_gout=0.99;
l=0.5;

 f1Y_B0=Y_B0 * (1-f_gout);

 selectHMX=1;
 Choice1=1;
 solverselect=1;
 Selector=2 ;
 
MAX = 100;
bcinterval=[-10^2 10^2];   
     
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Y_HMXprod=0.3;
Y_gTf=0;
Y_BF3 =0;
Y_CF2=0;
Y_C=0;
Y_C2=0;
Y_cTf=0;
 Y_B=0;
 fg=0.7;
 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;
  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, T_g)(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=@(T_g, Y_HMXprod, Y_gHMX)(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=@(T_g, Y_HMXprod, Y_gHMX)(Mg(T_g, Y_HMXprod, Y_gHMX).*P./(R.*T_g));
roc = ro_cTf*Y_cTf + ro_B*Y_B ;
ro_gf=@(T_g, Y_HMXprod, Y_gHMX)((1-fg)*roc+fg*rog(T_g, Y_HMXprod, Y_gHMX));

% теплоемкости газов______________________________________________________
c_NO2(T_g_<1200)=Pn5(16.108,75.89,-54.3874,14.30777,0.239423, T_g_(T_g_<1200))/M_NO2;
c_NO2(T_g_>=1200)=Pn5(56.82541,0.738053,-0.144721, 0.009777,-5.459911,T_g_(T_g_>=1200))/M_NO2;

c_NO(T_g_<1200)=Pn5(23.83, 12.589, -1.139, -1.497459, 0.214194, T_g_(T_g_<1200))/M_NO ;
c_NO(T_g_>=1200)=Pn5(35.9916912,0.957170,-0.148032, 0.009974,-3.004088, T_g_(T_g_>=1200))/M_NO  ;  

c_N2O(T_g_<1400)=(Pn5(27.67,51.148,-30.64, 6.847911, -0.157906, T_g_(T_g_<1400))/M_N2O );
c_N2O(T_g_>=1400)=Pn5(60.30274,1.034566,-0.192997, 0.012540, -6.860254,  T_g_(T_g_>=1400))/M_N2O  ;

c_H2O(T_g_<1700)=Pn5(30.092,6.832514,6.793435,-2.534480,0.082139,T_g_(T_g_<1700) )/M_H2O;
c_H2O(T_g_>=1700)=Pn5(41.96426,8.622053,-1.499780,0.098119,-11.15764, T_g_(T_g_>=1700) )/M_H2O;

c_CH2O(T_g_<1200)=Pn5(5.19,93.2,-44.85, 7.882279,0.551175, T_g_(T_g_<1200))/M_CH2O ;
c_CH2O(T_g_>=1200)=Pn5(71.35268,6.174497,-1.191090, 0.079564,-15.58917,T_g_(T_g_>=1200) )/M_CH2O ;

c_CO(T_g_<1300)=Pn5(25.56759,6.096130,4.054656, -2.671301, 0.131021, T_g_(T_g_<1300))/M_CO ;
c_CO(T_g_>=1300)=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) = Pn5(21.28631, 130.3006, -109.9919, 34.28838, -0.073386, T_g_(T_g_<1000))/ M_BF3;
C_BF3(T_g_>=1000) =  Pn5(81.23696, 1.096330, -0.226830, 0.015981, -6.366625, T_g_(T_g_>=1000))/M_BF3;

c_gTf(T_g_<1100)=Pn5(43.55126,175.9079,-138.7331,40.35619,-0.381260,  T_g_(T_g_<1100))/M_Tf;
c_gTf(T_g_>=1100)=Pn5(129.9776,1.690918,-0.340087,0.023448,-10.83204,  T_g_(T_g_>=1100))/M_Tf;

C_CF2(T_g_<600) =Pn5(8.825772, 125.3652, -129.4940, 50.33101, 0.259749, T_g_(T_g_<600))/ M_CF2;
C_CF2(T_g_>=600) =  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) = Pn5(123.9092,-348.0067,485.0971, -232.7994, -1.240298, T_g_(T_g_<700))/M_C2;
C_C2(T_g_>=700) = 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=@(T_g, Y_HMXprod, Y_gHMX)(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.04;
c_cTf(T_g_>=TmTf)= (0.61488+0.001194*1000*tau(T_g_));

c_B = @(T_g)(10.18574 + 29.24415*tau(T_g) -18.02137*tau(T_g).^2 +4.212326*tau(T_g).^3 )/10.8 ;    
cc=@(T_g)(c_cTf(T_g).*Y_cTf+ c_B(T_g).*Y_B); 
c_gf=@(T_g, Y_HMXprod, Y_gHMX)(((1-fg).*roc(T_g, Y_HMXprod, Y_gHMX).*cc(T_g)+fg.*rog(T_g, Y_HMXprod, Y_gHMX).*cg(T_g, Y_HMXprod, Y_gHMX))./ro_gf(T_g, Y_HMXprod, Y_gHMX));

% теплопроводности_________________________________________________________
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=@(T_g, Y_HMXprod, Y_gHMX)(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= @(T_g, Y_HMXprod, Y_gHMX)( ((1-fg).*lc.*cc(T_g)+fg.*lg(T_g, Y_HMXprod, Y_gHMX).*cg(T_g, Y_HMXprod, Y_gHMX))./((1-fg).*cc(T_g, Y_HMXprod, Y_gHMX)+fg.*cg(T_g, Y_HMXprod, Y_gHMX)));
 %________________________________________________________________________
% кинетмка
%Tf->C2F4
% C2F4->2CF2
% B+3CF2->2BF3+3C
% 2C->C2
switch selectHMX
    case 1
G_gHMX = @(T_g, Y_HMXprod, Y_gHMX)(fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_gHMX.*Ar(10^14.2, 165268, T_g)/M_HMX); %
    case 2
 G_gHMX =fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_gHMX.*Ar(10^12.5, 158992)/M_HMX;
    case 3
 G_gHMX =fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_gHMX.*Ar(5*10^13, 158992)/M_HMX;        
  
end
G_Tf= @(T_g, Y_HMXprod, Y_gHMX)(((1-fg).*ro_cTf.*Y_cTf.*Ar(kTf,ETf, T_g)/M_Tf));
G_CF2dec = @(T_g, Y_HMXprod, Y_gHMX)(fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_gTf.*T_g.^(-0.5).*Ar(7.82*10^15,2.33*10^5, T_g)/M_C2F4);
G_B_CF2 =  @(T_g, Y_HMXprod, Y_gHMX)(fg*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_CF2.*(Y_B*ro_B).^3.*T_g.^(0.5).*Ar(4*10^14,8.37*10^4, T_g)/(M_CF2*M_B^3));
G_C = @(T_g, Y_HMXprod, Y_gHMX)(fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).^2.*Y_C.^2.*T_g.^(-1.6).*Ar(1.80*10^21,0, T_g)/M_C);

w_gHMX=  @(T_g, Y_HMXprod, Y_gHMX)(-M_HMX*G_gHMX(T_g, Y_HMXprod, Y_gHMX));
w_gHMXprod = @(T_g, Y_HMXprod, Y_gHMX)(M_HMX*G_gHMX(T_g, Y_HMXprod, Y_gHMX));
%Tf->C2F4
w_cTf=@(T_g, Y_HMXprod, Y_gHMX)(-M_Tf*G_Tf);
 w_gTf =@(T_g, Y_HMXprod, Y_gHMX)( -M_Tf*w_cTf);
% B+3CF2->BF3+3C
w_B=@(T_g, Y_HMXprod, Y_gHMX)(-M_B*G_B_CF2);
% C2F4->2CF2, B+3CF2->2BF3+3C
w_CF2=@(T_g, Y_HMXprod, Y_gHMX)(2*M_CF2*G_CF2dec-3*M_CF2*G_B_CF2);
% B+3CF2->2BF3+3C
w_BF3 =@(T_g, Y_HMXprod, Y_gHMX)(2*M_BF3* G_B_CF2);
% B+3CF2->2BF3+3C, 2C->C2
w_C = @(T_g, Y_HMXprod, Y_gHMX)(-2*M_C*G_C +3*M_C*G_B_CF2); 
w_C2 =@(T_g, Y_HMXprod, Y_gHMX)( M_C2*G_C);
G_condphase= @(T_g, Y_HMXprod, Y_gHMX)(-w_cTf-w_B);
% тепловыделение___________________________________________________________
Q_react_g=@(T_g, Y_HMXprod, Y_gHMX)((Q_gHMX*w_gHMX(T_g, Y_gHMX)/(wm*dHm))); %0*( -w_cTf*824/M_Tf - w_gTf*658.56/M_Tf + 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=@(T_g, Y_HMXprod, Y_gHMX)(cfm*m.*l./l_gf(T_g, Y_HMXprod, Y_gHMX)); 
Pe_dk = @(T_g, Y_HMXprod, Y_gHMX)(m*l./(D*rog(T_g, Y_HMXprod, Y_gHMX))); 
  %________________________________________________________________________
  res_T_g_ = zeros(n,1);
  res_Y_gHMX_ = zeros(n,1);
  res_Y_HMXprod_=zeros(n,1);
  %________________________________________________________________________
  res_T_g_(1) = T_g_(1)-T_fout;
  for i=2:n-1
  res_T_g_(i) =  (Pe_gm./Pe_g(T_g_(i), Y_HMXprod_(i),  Y_gHMX_(i)))*(T_g_(i+1)-2*T_g_(i)+T_g_(i-1))/dx^2 + ...
      0*(1./Pe_dk(T_g_(i), Y_HMXprod_(i),  Y_gHMX_(i))).*(c_gf(T_g_(i), Y_HMXprod_(i),  Y_gHMX_(i))/cfm)*ro_gf(T_g_(i), Y_HMXprod_(i),  Y_gHMX_(i)).*((T_g_(i+1)-T_g_(i-1))/(2*dx)).*Y_gHMX_(i)*(Y_gHMX_(i+1)-Y_gHMX_(i-1))/(2*dx) - ...
  (Pe_gm)^0.5 * (c_g(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i))/cfm)*(T_g_(i+1)-T_g_(i))/(2*dx) -  Q_react_g(T_g_(i), Y_HMXprod(i), Y_gHMX_(i));   %  Q_gHMX*Y_gHMX_(i).*ro_gf.*Ar(10^14.2, 165268) 1*
  end
  res_T_g_(n) = ( T_g_(n)-T_g_(n-1))/dx ;
  %________________________________________________________________________
  res_Y_gHMX_(1) = Y_gHMX_(1)-Y_gHMXout;
  for i = 2:n-1

    res_Y_gHMX_(i) =(1./Pe_dk(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i)))*(Y_gHMX_(i+1)-2*Y_gHMX_(i)+Y_gHMX_(i-1))/dx^2-(Y_gHMX_(i+1)-Y_gHMX_(i-1))/(2*dx)+w_gHMX(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i))/wm; %w_gHMX  Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
  end
  res_Y_gHMX_(n) =( Y_gHMX_(n)-Y_gHMX_(n-1))/dx;
  %________________________________________________________________________
  res_Y_HMXprod_(1) = Y_HMXprod_(1)-Y_HMXprod0;
 for i = 2:n-1
     
 res_Y_HMXprod_(i) =(1./Pe_dk(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i)))*(Y_HMXprod_(i+1)-2*Y_HMXprod_(i)+Y_HMXprod_(i-1))/dx^2-(Y_HMXprod_(i+1)-Y_HMXprod_(i-1))/(2*dx)+w_gHMXprod(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i))/wm; %w_gHMX  Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
  end
 res_Y_HMXprod_(n) =( Y_HMXprod_(n)-Y_HMXprod_(n-1))/dx;
  
  res = [res_T_g_;res_Y_gHMX_; res_Y_HMXprod_;];
end
end

[Cris LaPierre]

cc is a variable that got incorrectly set up as an anonymous function. I think it should be

cc=c_cTf.*Y_cTf+ c_B(T_g_).*Y_B; 

but you need to confirm that.

[Gleb]

unfortunately there is an error "??? Error using ==> plus

Matrix dimensions must agree.

Error in ==> Gas_system_4v3>fun at 191

cc=(c_cTf.*Y_cTf+ c_B(T_g_).*Y_B);"

[Cris LaPierre]

Not when I run it.

Check the dimensions of your variables to identify the differences. Then go back to where each of these variables were created to find any errors.

If you are not using it yet, I suggest learning how to use the debugger.

[Gleb]

Thank you! I fixed it.

but got a problem with

c_gf=@(T_g, Y_HMXprod, Y_gHMX)(((1-fg).*roc(T_g, Y_HMXprod, Y_gHMX).*cc+fg.*rog(T_g, Y_HMXprod, Y_gHMX).*cg(T_g, Y_HMXprod, Y_gHMX))./ro_gf(T_g, Y_HMXprod, Y_gHMX));

"

Index in position 1 exceeds array bounds (must not exceed 1).

Error in Gas_system_4v4>@(T_g,Y_HMXprod,Y_gHMX)(((1-fg).*roc(T_g,Y_HMXprod,Y_gHMX).*cc+fg.*rog(T_g,Y_HMXprod,Y_gHMX).*cg(T_g,Y_HMXprod,Y_gHMX))./ro_gf(T_g,Y_HMXprod,Y_gHMX)) (line 162)

c_gf=@(T_g, Y_HMXprod, Y_gHMX)(((1-fg).*roc(T_g, Y_HMXprod, Y_gHMX).*cc+fg.*rog(T_g, Y_HMXprod, Y_gHMX).*cg(T_g, Y_HMXprod, Y_gHMX))./ro_gf(T_g, Y_HMXprod, Y_gHMX));

"

Never met such one.

[Cris LaPierre]

Maybe a new error, but same issue. I think c_gf should not be an anonymous function. I think you are actually trying to create a vector of values.

[Cris LaPierre]

Perhaps a better way to approach this is start with the first anonymous functions your create (lines 135-148).

dHm=cfm;
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, T_g)(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=@(T_g, Y_HMXprod, Y_gHMX)(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=@(T_g, Y_HMXprod, Y_gHMX)(Mg(T_g, Y_HMXprod, Y_gHMX).*P./(R.*T_g));
roc = ro_cTf*Y_cTf + ro_B*Y_B ;
ro_gf=@(T_g, Y_HMXprod, Y_gHMX)((1-fg)*roc+fg*rog(T_g, Y_HMXprod, Y_gHMX));

We've already talked about tau, Pn2 and Pn5. They should remain anonymous functions. Take a hard look at the rest. Which do you think should be vectors of values instead of anonymous functions?

A good criteria to consider is if you ever create the inputs used when you call the function, or if you call the function multiple times with different input values.

[Gleb]

I think all functions in the for loop (Pe_dk, Pe_g, Q_react_g, w_gHMX) should be anonymos functions, because they change from i to i+1 (or from one space coordinate to another, because its a balance equations). I fixed most of the problems, but the error occurs

'Unable to perform assignment because the left and right sides have a different number of elements.

Error in Gas_system_4v4/fun (line 215)

res_T_g_(i) = (Pe_gm./Pe_g(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i)))*(T_g_(i+1)-2*T_g_(i)+T_g_(i-1))/dx^2 - ...'

I replaced Pe_g(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i))) and c_gf by dummy random constants and found out that code works and integration completes, So this functions spoils the code.

debugging shows that c_gf is n to n array, the same is about Pe_gm. So i think the problem lies here.

Here is my latest version

function Gas_system_4v4
clc
close all
n =50;
global T_fout Y_gHMXout Pe_g Pe_dk w_gHMX Pe_gm;
  x0 = [800*ones(n,1);1*ones(n,1);1*ones(n,1)];
  x = ((1:n)-1)/(n-1);
  Pe_gm=4.8157e+004;
 AlgoSelector=1;
 switch AlgoSelector
     case 1
 Algochise='levenberg-marquardt'; %  'trust-region-dogleg'
     case 2
 Algochise='trust-region-dogleg' ;
     case 3
 Algochise='trust-region' ;
 end
%  options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-5, 'FiniteDifferenceStepSize', 10^-6);
 % options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 5*10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'StepTolerance', 10^-7,  'FunctionTolerance', 1*10^-3, 'FiniteDifferenceStepSize', 10^-4);
   % 'OptimalityTolerance', 10^-11,
  %options=optimoptions('fsolve', 'Display','iter', 'MaxFunctionEvaluations', 10^6, 'MaxIterations', 10^6, 'Algorithm', Algochise, 'FunctionTolerance', 10^-4, 'OptimalityTolerance', 10^-3, 'StepTolerance', 10^-5  );    % optimoptions
  options=optimset( 'Display','iter', 'MaxFunEvals', 10^8, 'MaxIter', 10^6, 'Algorithm', Algochise, 'TolFun', 10^-7);   
  [sol,fval] = fsolve(@(x)fun(x,n),x0,  options);
  disp(fval);
  norm(fun(sol,n))
   figure
  plot(x*(  Pe_gm)^-0.5,sol(1:n))
  hold on
  figure
    plot(x,sol(n+1:2*n))
    hold on
     figure
    plot(x,sol(2*n+1:3*n))
    hold on
disp(Pe_g);
disp(Pe_dk);
function res = fun(Yg,n)
  x = ((1:n)-1)/(n-1);
  dx = 1/(n-1);
   T_fout =750;
Y_gHMXout =0.8;
 cfm=5*10^4;
  T_g_ = Yg(1:n);
  Y_gHMX_ = Yg(n+1:2*n);
  Y_HMXprod_ = Yg(2*n+1:3*n);
%   for i=1:length(T_g_)
% T_g= T_g_(i);
%   end
%    %for i=1:n
%   for i=1:length(Y_gHMX_) 
% Y_gHMX= Y_gHMX_(i);
%   end
%global T_fout  L_g Y_gHMXout  Y_B0 x  Q_react_g
% format long e
Y_HMXprod0=0.2;
f1Y_cTfout=0.0;
Y_B0 = 0.0;
 f_gout=0.99;
l=0.5;
 f1Y_B0=Y_B0 * (1-f_gout);
 selectHMX=1;
 Choice1=1;
 solverselect=1;
 Selector=2 ;
 
MAX = 100;
bcinterval=[-10^2 10^2];   
     
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Y_HMXprod=0.3;
Y_gTf=0;
Y_BF3 =0;
Y_CF2=0;
Y_C=0;
Y_C2=0;
Y_cTf=0;
 Y_B=0;
 fg=0.7;
 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;
  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, T_g)(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=@(T_g, Y_HMXprod, Y_gHMX)(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=@(T_g, Y_HMXprod, Y_gHMX)(Mg(T_g, Y_HMXprod, Y_gHMX).*P./(R.*T_g));
roc = ro_cTf*Y_cTf + ro_B*Y_B ;
ro_gf=@(T_g, Y_HMXprod, Y_gHMX)((1-fg)*roc+fg*rog(T_g, Y_HMXprod, Y_gHMX));
% теплоемкости газов______________________________________________________
c_NO2(T_g_<1200)=Pn5(16.108,75.89,-54.3874,14.30777,0.239423, T_g_(T_g_<1200))/M_NO2;
c_NO2(T_g_>=1200)=Pn5(56.82541,0.738053,-0.144721, 0.009777,-5.459911,T_g_(T_g_>=1200))/M_NO2;
c_NO(T_g_<1200)=Pn5(23.83, 12.589, -1.139, -1.497459, 0.214194, T_g_(T_g_<1200))/M_NO ;
c_NO(T_g_>=1200)=Pn5(35.9916912,0.957170,-0.148032, 0.009974,-3.004088, T_g_(T_g_>=1200))/M_NO  ;  
c_N2O(T_g_<1400)=(Pn5(27.67,51.148,-30.64, 6.847911, -0.157906, T_g_(T_g_<1400))/M_N2O );
c_N2O(T_g_>=1400)=Pn5(60.30274,1.034566,-0.192997, 0.012540, -6.860254,  T_g_(T_g_>=1400))/M_N2O  ;
c_H2O(T_g_<1700)=Pn5(30.092,6.832514,6.793435,-2.534480,0.082139,T_g_(T_g_<1700) )/M_H2O;
c_H2O(T_g_>=1700)=Pn5(41.96426,8.622053,-1.499780,0.098119,-11.15764, T_g_(T_g_>=1700) )/M_H2O;
c_CH2O(T_g_<1200)=Pn5(5.19,93.2,-44.85, 7.882279,0.551175, T_g_(T_g_<1200))/M_CH2O ;
c_CH2O(T_g_>=1200)=Pn5(71.35268,6.174497,-1.191090, 0.079564,-15.58917,T_g_(T_g_>=1200) )/M_CH2O ;
c_CO(T_g_<1300)=Pn5(25.56759,6.096130,4.054656, -2.671301, 0.131021, T_g_(T_g_<1300))/M_CO ;
c_CO(T_g_>=1300)=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) = Pn5(21.28631, 130.3006, -109.9919, 34.28838, -0.073386, T_g_(T_g_<1000))/ M_BF3;
C_BF3(T_g_>=1000) =  Pn5(81.23696, 1.096330, -0.226830, 0.015981, -6.366625, T_g_(T_g_>=1000))/M_BF3;
c_gTf(T_g_<1100)=Pn5(43.55126,175.9079,-138.7331,40.35619,-0.381260,  T_g_(T_g_<1100))/M_Tf;
c_gTf(T_g_>=1100)=Pn5(129.9776,1.690918,-0.340087,0.023448,-10.83204,  T_g_(T_g_>=1100))/M_Tf;
C_CF2(T_g_<600) =Pn5(8.825772, 125.3652, -129.4940, 50.33101, 0.259749, T_g_(T_g_<600))/ M_CF2;
C_CF2(T_g_>=600) =  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) = Pn5(123.9092,-348.0067,485.0971, -232.7994, -1.240298, T_g_(T_g_<700))/M_C2;
C_C2(T_g_>=700) = 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=@(T_g, Y_HMXprod, Y_gHMX)(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.04;
c_cTf(T_g_>=TmTf)= (0.61488+0.001194*1000*tau(T_g_));
c_B = @(T_g)(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(T_g_).*Y_B; 
c_gf=(((1-fg).*roc.*cc+fg.*rog(T_g_, Y_HMXprod_, Y_gHMX_).*cg(T_g_, Y_HMXprod_, Y_gHMX_))./ro_gf(T_g_, Y_HMXprod_, Y_gHMX_));

% теплопроводности_________________________________________________________
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=@(T_g, Y_HMXprod, Y_gHMX)(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= @(T_g, Y_HMXprod, Y_gHMX)( ((1-fg).*lc.*cc+fg.*lg(T_g, Y_HMXprod, Y_gHMX).*cg(T_g, Y_HMXprod, Y_gHMX))./((1-fg).*cc+fg.*cg(T_g, Y_HMXprod, Y_gHMX)));
 %________________________________________________________________________
% кинетмка
%Tf->C2F4
% C2F4->2CF2
% B+3CF2->2BF3+3C
% 2C->C2
switch selectHMX
    case 1
G_gHMX = @(T_g, Y_HMXprod, Y_gHMX)(fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_gHMX.*Ar(10^14.2, 165268, T_g)/M_HMX); %
    case 2
 G_gHMX =fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_gHMX.*Ar(10^12.5, 158992)/M_HMX;
    case 3
 G_gHMX =fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_gHMX.*Ar(5*10^13, 158992)/M_HMX;        
  
end
G_Tf= @(T_g, Y_HMXprod, Y_gHMX)(((1-fg).*ro_cTf.*Y_cTf.*Ar(kTf,ETf, T_g)/M_Tf));
G_CF2dec = @(T_g, Y_HMXprod, Y_gHMX)(fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_gTf.*T_g.^(-0.5).*Ar(7.82*10^15,2.33*10^5, T_g)/M_C2F4);
G_B_CF2 =  @(T_g, Y_HMXprod, Y_gHMX)(fg*ro_gf(T_g, Y_HMXprod, Y_gHMX).*Y_CF2.*(Y_B*ro_B).^3.*T_g.^(0.5).*Ar(4*10^14,8.37*10^4, T_g)/(M_CF2*M_B^3));
G_C = @(T_g, Y_HMXprod, Y_gHMX)(fg.*ro_gf(T_g, Y_HMXprod, Y_gHMX).^2.*Y_C.^2.*T_g.^(-1.6).*Ar(1.80*10^21,0, T_g)/M_C);
w_gHMX=  @(T_g, Y_HMXprod, Y_gHMX)(-M_HMX*G_gHMX(T_g, Y_HMXprod, Y_gHMX));
w_gHMXprod = @(T_g, Y_HMXprod, Y_gHMX)(M_HMX*G_gHMX(T_g, Y_HMXprod, Y_gHMX));
%Tf->C2F4
w_cTf=@(T_g, Y_HMXprod, Y_gHMX)(-M_Tf*G_Tf);
 w_gTf =@(T_g, Y_HMXprod, Y_gHMX)( -M_Tf*w_cTf);
% B+3CF2->BF3+3C
w_B=@(T_g, Y_HMXprod, Y_gHMX)(-M_B*G_B_CF2);
% C2F4->2CF2, B+3CF2->2BF3+3C
w_CF2=@(T_g, Y_HMXprod, Y_gHMX)(2*M_CF2*G_CF2dec-3*M_CF2*G_B_CF2);
% B+3CF2->2BF3+3C
w_BF3 =@(T_g, Y_HMXprod, Y_gHMX)(2*M_BF3* G_B_CF2);
% B+3CF2->2BF3+3C, 2C->C2
w_C = @(T_g, Y_HMXprod, Y_gHMX)(-2*M_C*G_C +3*M_C*G_B_CF2); 
w_C2 =@(T_g, Y_HMXprod, Y_gHMX)( M_C2*G_C);
G_condphase= @(T_g, Y_HMXprod, Y_gHMX)(-w_cTf-w_B);
% тепловыделение___________________________________________________________
Q_react_g=@(T_g, Y_HMXprod, Y_gHMX)((Q_gHMX*w_gHMX(T_g, Y_HMXprod, Y_gHMX)/(wm*dHm))); %0*( -w_cTf*824/M_Tf - w_gTf*658.56/M_Tf + 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=@(T_g, Y_HMXprod, Y_gHMX)(cfm*m.*l./l_gf(T_g, Y_HMXprod, Y_gHMX)); 
Pe_dk = @(T_g, Y_HMXprod, Y_gHMX)(m*l./(D*rog(T_g, Y_HMXprod, Y_gHMX))); 
  %________________________________________________________________________
  res_T_g_ = zeros(n,1);
  res_Y_gHMX_ = zeros(n,1);
  res_Y_HMXprod_=zeros(n,1);
  %________________________________________________________________________
  res_T_g_(1) = T_g_(1)-T_fout;
  for i=2:n-1
 
  res_T_g_(i) =  (Pe_gm./Pe_g(T_g_(i), Y_HMXprod_(i),  Y_gHMX_(i)))*(T_g_(i+1)-2*T_g_(i)+T_g_(i-1))/dx^2 - ...
 (Pe_gm)^0.5 * (c_gf/cfm)*(T_g_(i+1)-T_g_(i))/(2*dx) -  Q_react_g(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i));   %  Q_gHMX*Y_gHMX_(i).*ro_gf.*Ar(10^14.2, 165268) 1* 0*(1./Pe_dk(T_g_(i), Y_HMXprod_(i),  Y_gHMX_(i))).*(c_gf(T_g_(i), Y_HMXprod_(i),  Y_gHMX_(i))/cfm)*ro_gf(T_g_(i), Y_HMXprod_(i),  Y_gHMX_(i)).*((T_g_(i+1)-T_g_(i-1))/(2*dx)).*Y_gHMX_(i)*(Y_gHMX_(i+1)-Y_gHMX_(i-1))/(2*dx) 
 end
  res_T_g_(n) = ( T_g_(n)-T_g_(n-1))/dx ;
  %________________________________________________________________________
  res_Y_gHMX_(1) = Y_gHMX_(1)-Y_gHMXout;
  for i = 2:n-1
    res_Y_gHMX_(i) =(1./Pe_dk(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i)))*(Y_gHMX_(i+1)-2*Y_gHMX_(i)+Y_gHMX_(i-1))/dx^2-(Y_gHMX_(i+1)-Y_gHMX_(i-1))/(2*dx)+w_gHMX(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i))/wm; %w_gHMX  Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
  end
  res_Y_gHMX_(n) =( Y_gHMX_(n)-Y_gHMX_(n-1))/dx;
  %________________________________________________________________________
  res_Y_HMXprod_(1) = Y_HMXprod_(1)-Y_HMXprod0;
 for i = 2:n-1
     
 res_Y_HMXprod_(i) =(1./Pe_dk(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i)))*(Y_HMXprod_(i+1)-2*Y_HMXprod_(i)+Y_HMXprod_(i-1))/dx^2-(Y_HMXprod_(i+1)-Y_HMXprod_(i-1))/(2*dx)+w_gHMXprod(T_g_(i), Y_HMXprod_(i), Y_gHMX_(i))/wm; %w_gHMX  Y_gHMX_(i).*ro_gf.*(10^14.2)*exp(-165268./(8.314*T_g_(i) ))
  end
 res_Y_HMXprod_(n) =( Y_HMXprod_(n)-Y_HMXprod_(n-1))/dx;
  
  res = [res_T_g_;res_Y_gHMX_; res_Y_HMXprod_;];
end
end

[Cris LaPierre]

Let's focus on the ones in the code I shared for now.

As for the for loops at the bottom, you don't need them. I would use indexing.

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);

I couldn't help but notice that a large section of your calculation of res_T_g_ is multiplied by 0. Since this doesn't contribute anything to the calculation, you could simply the code by removing it.

res_T_g_(ind) =  (Pe_gm./Pe_g(ind)).*(T_g_(ind+1)-2*T_g_(ind)+T_g_(ind-1))/dx^2 - ...
    (Pe_gm)^0.5 .* (c_gf(ind)/cfm).*(T_g_(ind+1)-T_g_(ind))/(2*dx) -  Q_react_g(ind);

[Gleb]

Well yes, this section can be removed, for now im solving the simplest model.

indexing is an intriguing move, is there some advantage over traditional for loop?""

"As for the for loops at the bottom, you don't need them". Well... There are 3 equations, and 3 loops, isnt it?

[Cris LaPierre]

I've shown you how to remove all the for loops, so I think you only need 4 anonymous functions in your code.

 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_)));

See if you can rewrite your code to remove all other anonymous functions.

[Cris LaPierre]

is there some advantage over traditional for loop?

Loops = slower code. In MATLAB, you can often avoid them by taking advantage of elementwise operators.

Well... There are 3 equations, and 3 loops, isnt it?

I've got to leave something for you to do. You can use the example I shared to see how to remove the loops from the remaining 2 equations.