Index exceeds matrix dimensions.

clearvars
close all
clc

global eta eps rext rint delta rhoc cpi PinSweep
global DH0 Trif mui ki MW S R Ac Dp p ls Tw

% Gas Constant [J/(kmol*K)]

R = 8314.45;

% Catalyst Density [kgcat/m^3]

rhoc = 2400;

% Membrane Thickness [m]

delta = 5e-6;

% Reaction efficiency [-]

eta = 0.28;

%External Radius [m]

rext = 0.08; 

% Internal Radius [m]

rint = 0.045;

% Cross Sectional Area [m^2]

Ac = pi*rint^2;

% Reactor Length [m]

L = 10;

% Void Fraction [-]

eps = 0.6;

% Stoichiometric Matrix [-]

nu = [-1    -1      1      1      0      0      0];
      %CO   %H2O   %CO2    %H2  %Inert %H2perm %sweep

% Stoichiometric Matrix at Permeate side [-]

nup = [0      0      0      -1      0     -1      0];
      %CO   %H2O   %CO2    %H2  %Inert %H2perm %sweep

% Specific Heats matrix [kJ/kmol*K]

cpi = 1e-3*[2.9108e+04   3.3363e+04   2.9370e+04   2.7617e+04   2.9105e+04    2.7617e+04   3.3363e+04     %C1
            8.7730e+03   2.6790e+04   3.4540e+04   9.5600e+04   8.6149e+03    9.5600e+04   2.6790e+04     %C2
            3.0851e+03   2.6105e+03   1.4280e+03   2.4660e+03   1.7016e+03    2.4660e+03   2.6105e+03     %C3
            8.4553e+03   8.8960e+03   2.6400e+04   3.7600e+03   1.0347e+02    3.7600e+03   8.8960e+03     %C4
            1.5382e+03   1.1690e+03   5.8800e+02   5.6760e+02   9.0979e+02    5.6760e+02   1.1690e+03];   %C5
               %CO          %H2O         %CO2         %H2          %Inert       %H2perm       %sweep  

% Viscosity [Pa*s]

mui = [ 1.1127e-06    1.7096e-08      2.148e-06      1.797e-07   6.5592e-07       1.797e-07     1.7096e-08    %C1
        0.5338        1.1146          0.46           0.685       0.6081           0.685         1.1146        %C2
       94.7           0             290             -0.59        54.714          -0.59          0             %C3
        0             0               0            140           0              140             0        ];   %C4

            %CO            %H2O         %CO2         %H2           %Inert         %H2perm         %sweep 

% Thermal conductivity [W/m*K]

ki = [  5.9882e-04     6.2041e-06     3.69         2.653e-03  3.3143e-04      2.653e-03     6.2041e-06     %C1
        0.6863         1.3973        -0.3838       0.7452     0.7722          0.745         1.3973         %C2
       57.13           0            964           12          16.323         12             0              %C3
      501.92           0              1.86e+06     0          373.72          0             0          ];  %C4

            %CO         %H2O         %CO2           %H2         %Inert        %H2perm       %sweep 

% Ideal gas enthalpy of formation [kJ/kmol]  Trif =298.15 K

DH0 = 1e-3*[-1.1053e+08   -2.41814e+08    -3.9351e+08        0            0            0       -2.41814e+08 ];
                %CO            %H2O           %CO2          %H2        %Inert       %H2perm       %sweep 

% Reference Temperature [K]

Trif = 298.15;

% Wall Temperature [K]

Tw = 628.15;

% Initial Sweep Temperature [K]

TinSweep = 628.15;

% Initial Temperature [K]

Tin = 628;

% Initial Sweep Pressure [kPa]

PinSweep = 101;

% Initial Pressure [kPa]

Pin = 2376.53;

% Initial Sweep Flow Rate [kmol/h] 

FinSweep = 8;

% Initial Flow Rate [kmol/h]

Ftot0 = 14;

% Inlet Compositions at Permeation Side [-]

F0Perm = 0;

% Inlet Compositions of reacting Mixture [-]

x0 = [0.0521   0.4453   0.0272   0.4649   0.0105];
        %CO     %H2O     %CO2      %H2    %Inert

% MW [kg/kmol]

MW = [ 28.010        18.015       44.010      2.016       28.013       2.016      18.015 ];
         %CO          %H2O         %CO2        %H2        %Inert      %H2perm      %sweep 

% Normal Boiling Temperature and corresponding Parameters [K]

Tb = [81.15   373.15   194.65   79    77.35      79      373.15 ];
       %CO     %H2O     %CO2    %H2   %Inert   %H2perm   %sweep       
S = 1.5.*Tb; 

% Global Heat Exchange coefficient for permeate side [kJ/(m^2*h*K)]

Um = 8.64;

% Particle Diameter [m]

Dp = 0.0053; 

% Thermal Conductivity of packing Material [kJ/(m*h*K)]

ls = 1.256;

% Solid Surface Emissivity [-]

p = 0.8;

% Inlet Conditions 

In = [Ftot0*x0, F0Perm, FinSweep, Pin, Tin, TinSweep - 100]; 

% Sistem Solution

zspan = [0 L]; % [m]
options = odeset('Reltol',1e-11,'Abstol',1e-12);
[Z,X] = ode15s(@Tronky,zspan,In,options);

function f = Tronky(z,x)

% Variable Definition

global eta eps rext rint delta rhoc cpi PinSweep
global DH0 Trif mui ki MW S R Ac Dp p ls Tw

Fi = x(1:5);
FiPS = x(6:7);
P = x(8);
T = x(9);
Tperm = x(10);

% Other Variables

Ftot = sum(Fi); % [kmol/h]
xi = Fi./Ftot; % [-]
Pi = P*xi; % [kPa]
FtotPS = sum(FiPS); % [kmol/h]
xiPS = FiPS./FtotPS; % [-]
PiPS = PinSweep*xiPS; % [kPa]
Tmem = (T + Tperm)/2; % [K]

% Reaction, Adsorption, Equilibrium Constants

k = exp(-29364/(1.987*T) + 40.32/1.987); % [kmol/(kgcat*h)]
Kco = exp(3064/(1.987*T) - 6.74/1.987)*1/(101325); % [1/Pa]
Kh2o = exp(-6216/(1.987*T) + 12.77/1.987); % [1/Pa]
Kco2 = exp(12542/(1.987*T) - 18.45/1.987); % [1/Pa]
Keq = exp(4400/T - 4.063); % [-]

% Reaction Rate [kmol/m^3*h]

rco = k*Kco*Kh2o*(Pi(1)*Pi(2) - Pi(3)*Pi(4)/Keq)*(1 + Kco*Pi(1) + Kh2o*Pi(2) + Kco2*Pi(4))^-2*rhoc; 
Rco = rco*eta*(1 - eps)*pi*(rext^2 - rint^2); 

% Membrane Permeability [kmol/(m*h*Pa^0.5)]

Bh = 2.95*1e-4*3600/1e3*exp(-5833.5/Tmem);

% Hydrogen Flux [kmol/(h*m^2)]

Jh2perm = Bh/delta*(Pi(4)^0.5 - PiPS(2)^0.5);
Rh2 = Jh2perm*2*pi*rint;

% Specific Heats, Enthalpies and molecular Weights of reacting Mixture

for i = 1:5
    cp(i) = cpi(1,i) + cpi(2,i)*((cpi(3,i)/T)/sinh(cpi(3,i)/T))^2 +...
        cpi(4,i)*((cpi(5,i)/T)/cosh(cpi(5,i)/T))^2; % [kJ/(kmol*K)]
    DHi(i) = DH0(i) + (2*cpi(2,i)*cpi(3,i))*((1/(exp(2*cpi(3,i)/T) - 1) - 1/(exp(2*cpi(3,i)/Trif) - 1)) + ...
        (1/(exp(2*cpi(3,i)/T) + 1) - 1/(exp(2*cpi(3,i)/Trif) + 1))); % [kJ/kmol]
    cpmixi(i) = cp(i)*xi(i); % [kJ/(mol*K)]
    MWi(i) = MW(i)*xi(i);
end

% Mixture Molecular Weight [kg/kmol]

MWm = sum(MWi);

% Mixture Density [kg/m^3]

rhom = P/(R*T)*MWm;

% Mixture Specific Heat [kJ/(mol*K)]

cpm = sum(cpmixi(i));

% Reaction Heat [kJ/mol]

DHr = -DHi(1) - DHi(2) + DHi(3) + DHi(4);

% Viscosities and Conductivities 

for i = 1:5

      mu = mui(1,i)*(T^(mui(2,i)))/(1 + mui(3,i)/T + mui(4,i)/(T^2)); % [Pa*s]
      k = 3600/1e3*(ki(1,i)*(T^(ki(2,i)))/(1 + ki(3,i)/T + ki(4,i)/(T^2))); % [kJ/(m*h*K)]
  end

% Mixing Rules for Viscosities and Conductivities

for i = 1:5
    for j= 1:5

          phi(i,j) = (1 + ((mu(i)./mu(j))^(1/2)).*((MW(j)./MW(i))^(1/4))).^2./...
              (8.*(1 + MW(i)./MW(j))).^(1/2); % [-]
          Ss(i,j) = 0.735*sqrt(S(i)*S(j)); % [K]
          Ai(i,j) = 0.25*(1 + (mu(i)/mu(j)*(MW(j)/MW(i))^(3/4)*((T + S(i))/(T + S(j))))^(1/2))^2*...
              ((T + Ss(i,j))/(T + S(i))); % [-]
          DENm(i) = phi(i,j)*xi(j); % [-]
          DENk(i) = Ai(i,j)*xi(j); % [-]
      end
      DENmt(i) = sum(DENm); % [-]
      DENkt(i) = sum(DENk); % 
      NOMm(i) = mu(i)*xi(i); % [Pa*s]
      NOMk(i) = k(i)*xi(i); % [kJ/(m*h*K)]

end
mum = sum(NOMm)/sum(DENmt);
km = sum(NOMk)/sum(DENkt);

% Mixture Velocity [m/s]

v = Ftot*MWm*1/(rhom*Ac*3600);

% Adimensional Numbers [-]
Re = rhom*v*Dp/mum;
Pr = mum*cpm/km;

% Parameters for global Heat Exchange Coefficient

alpharu = (0.8171*(T/100)^3)/(1 + eps/(2*(1 - eps))*((1 - p)/p));
alphars = 0.8171*(p/(2 - p))*(T/100)^3;
ler0 = eps*(km + 0.95*alpharu*Dp) + ((0.95*(1 - eps))/(2/(3*ls) + 1/(10*km + alphars*Dp))); % [kJ/(m*K*h)]
ler = ler0 + 0.111*km*(Re*Pr^(1/3))/(1 + 46*(Dp/(2*rext))^2); % [kJ/(m*K*h)]
alphaw = 8.694*(ler0/(2*rext)^(4/3)) + 0.512*km*2*rext*Re*Pr^(1/3)/Dp;
Bi = alphaw*rext/ler; % [-]
U = (1/alphaw + rext/(3*ler)*(Bi + 3)/(Bi + 4))^(-1); % [kJ/(m^2*K*h)]

% Specific Heats and Enthalpies of Permeation Side

for i = 6:7

      cp(i) = cpi(1,i) + cpi(2,i)*((cpi(3,i)/T)/sinh(cpi(3,i)/T))^2 +...
          cpi(4,i)*((cpi(5,i)/T)/cosh(cpi(5,i)/T))^2; % [kJ/(kmol*K)]
      DHi(i) = DH0(i) + (2*cpi(2,i)*cpi(3,i))*((1/(exp(2*cpi(3,i)/T) - 1) - 1/(exp(2*cpi(3,i)/Trif) - 1)) + ...
          (1/(exp(2*cpi(3,i)/T) + 1) - 1/(exp(2*cpi(3,i)/Trif) + 1))); % [kJ/kmol]
      DHTpermi(i) = DH0(i) + (2*cpi(2,i)*cpi(3,i))*((1/(exp(2*cpi(3,i)/Tperm) - 1) - 1/(exp(2*cpi(3,i)/Trif) - 1)) + ...
          (1/(exp(2*cpi(3,i)/Tperm) + 1) - 1/(exp(2*cpi(3,i)/Trif) + 1))); % [kJ/kmol]
  end

% Ergun Expression

fv = 150 + 1.75*Re/(1 - eps);
Pressure = fv*v*mum/Dp^2*((1 - eps)^2/eps^3);

% Balance Equations

f(1) = -Rco;
f(2) = -Rco;
f(3) = Rco;
f(4) = Rco - Rh2;
f(5) = 0;
f(6) = -Rh2;
f(7) = 0;
f(8) = Pressure;
f(9) = 1/(Ftot*cpm)*(Rco*DHr + U*2*pi*rext*(Tw - T) - Um*2*pi*rint*(T - Tperm));
f(10) = 1/(FiPS(1)*cp(6) + FiSweep*cp(7))*(Um*2*pi*rint*(T - Tperm) + Rh2*(DHTpermi(6) - DHi(6)));

f = f(:);

end

Answer 1

[Rik]

Let's start with the obvious problem: You didn't post a question, just a wall of unformatted code. Ask an actual question instead of letting us guess what you want us to do.

Another issue is the use of global variables. Global variables are a bad idea. In some very specific cases they are unavoidable, but under almost every circumstance you will be able to solve your problem without them. Use function arguments to pass variables from one workspace to another. If you really, absolutely, undoubtedly need global variables, use as few as possible, and with long names that extensively describe what they are and what they should do. This minimizes the risk of using the same global variable name in two scripts/functions.

After spending some time pasting your code into Matlab and adding back some returns, I have many m-lint warnings. Sometimes they can be ignored, but it is always a good place to start if you encounter problems. No warnings is no guarantee for good code, but at least you have a solid foundation.

This looks like a lot of manual labour, where the chance of a typo is huge. Finding out which index is wrong might be a needle in a haystack. In this case however, the solution is with my previous remark: there is an m-lint warning about S not being unused or unset. This means that it is an empty matrix when you try to access elements 1 through 5.

Next time, have a read here and here. It will greatly improve your chances of getting an answer.