Index exceeds array bounds in BVP problem

I'm trying to make FEM for my boundary value problem. matlab gives an error "Index exceeds array bounds.

Error in Gas_system_4>fun (line 28)

Y_gHMX= Y_gHMX_(i);

"

I assume its because i use element of vector Y_gHMX_ outside of the loop. But how to fix it

function Gas_system_4
clc
close all
n = 100;
  x0 = ones(2*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(Yg,n)
%  some initial and boundary conditions and constanta
    % ...
    % ...
  x = ((1:n)-1)/(n-1);
  dx = 1/(n-1);

  T_g_ = Yg(1:n);
  Y_gHMX_ = Yg(n+1:2*n);
  for i=1:n
T_g= T_g_(i);
  end
  for i=n+1:2*n
Y_gHMX= Y_gHMX_(i);
  end
  %  some functions of T_g and  Y_gHMX
 Pe_g=... 
Pe_dk = ..
  w_gHMX=...
  Q_react_g=...
  
  
  res_T_g_ = zeros(n,1);
  res_Y_gHMX_ = zeros(n,1);

  
  res_T_g_(1) = T_g_(1)-T_fout;
  for i=2:n-1
    res_T_g_(i) = (1./Pe_g).*(T_g_(i+1)-2*T_g_(i)+T_g_(i-1))/dx^2 - (T_g_(i+1)-T_g_(i-1))/(2*dx) + Q_react_g;
  end
  res_T_g(n) = T_g_(n);
  res_Y_gHMX_(1) = Y_gHMXout(1);
  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;
  end
  res_Y_gHMX(n) = Y_gHMX_(n)-0;

  res = [res_T_g;res_Y_gHMX]; 
end

Answer 1

[Cris LaPierre]

Because i exceeds the length of Y_gHMX.

I see your index relates to the values you extracted from Y_gHMX_. However, once extracted, the index of Y_gHMX starts at 1, not n+1, and ends at (2*n-(n+1)+1). You could just use length(Y_gHMX) instead.

[Gleb]

Thank you!