Matlab - Problema de valor Inicial ODE113

Como estan comunidad,

hace un par de dias ando batallando con este codigo para resolver ecuaciones diferenciales ordinarias. El problema se trata de un balance de masa dinamico en funcion del tiempo: y a partir de x1 y x2 estas se consumiran y generaran como productos x3 hasta x12

lo que pasa es q este problema especificando

% Condiciones iniciales
x0=[0.25,0.75,0,0,0,0,0,0,0,0,0,0]; %Vector inicial

solo me genera resultados de x3 y x4.

la formacion de x5 hasta x12, depende mas que todo de la formacion de x3 ...si bien mi problema me genera de valores de x3--> esta cantidad no se transforma en x5 hasta x12

Pero cuando, cambio las condiciones iniciales a un vector inicial:
x0=[0.25,0.75,0 .075,0 .075,0,0,0,0,0,0,0,0]; %Vector inicial

si se presenta una generacion de resultados para las variables x5 hasta x12.

Como puede pasar esto?, si es que mi problema esta diseniado para generar x3 desde que se empieza a resolver el problema. Alguien podria ayudarme porfavor?

Aqui les envio el Script de matlab

Muchas gracias de antemano

__________________________________________________________________________________________

clear
clc
close all
%Parametros
R=8.3144598; %j/mol*°K

% Parametros
T=320+273.15; %Temperatura en °K
pt=0.2; % Presion im MPa
WHSV=1; %h^-1

% Condiciones iniciales
x0=[0.25,0.50,0.075,0.075,0,0,0,0,0,0,0,0]; %Vector inicial
x=x0;
tf=1/WHSV; %Zeit am ende in h

%Masa Molar %g/mol
M1=44.0095; %CO2
M2=2.01588; %H2
M3=28.0101; %CO
M4=18.01528; %H2O
M5=16.04246; %CH4
M6=28.05316; %C2H4
M7=42.07974; %C3H6
M8=44.09562; %C3H8
M9=56.10632; %C4H8
M10=70.1329; %C5H10
M11=140.2658; %C10H20
M12=88.10512; %C4H8O2

% syms x Masa molar total de la mezcla en funcion de las variables resultado (Fracciones molares)
%MT=((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))

%Constantes de velocidad de reaccion a temperatura T
k1=@(t)(5.5E6)*exp(-72.2E3/(R*T)); % mol/(g*h*MPa^2)
k2=@(t)(1.1E9)*exp(-98E3/(R*T)); % mol/(g*h*MPa^1.2)
k3=@(t)(1.2E4)*exp(-45E3/(R*T)); % mol/(g*h*MPa^1.2)
k4=@(t)(1.2E5)*exp(-53.6E3/(R*T)); % mol/(g*h*MPa^1.2)
k5=@(t)(1.8E4)*exp(-46.4E3/(R*T)); % mol/(g*h*MPa^1.2)
k6=@(t)(1.1E3)*exp(-31.5E3/(R*T)); % mol/(g*h*MPa^1.2)
k7=@(t)(9.0E2)*exp(-31.5E3/(R*T)); % mol/(g*h*MPa^1.2)
k8=@(t)(4E3)*exp(-33.4E3/(R*T)); % mol/(g*h*MPa^1.2)
k9=@(t)(3.6E3)*exp(-31.5E3/(R*T)); % mol/(g*h*MPa^1.2)
k10=@(t)(2.7E3)*exp(-25.8E3/(R*T)); % mol/(g*h*MPa^1.3)
k11=@(t)(1.4E4)*exp(-35.5E3/(R*T)); % mol/(g*h*MPa^1.8)
k12=@(t)(3.7E4)*exp(-40.7E3/(R*T)); % mol/(g*h*MPa^1.6)
k13=@(t)(6.0E7)*exp(-88.8E3/(R*T)); % mol/(g*h*MPa^1.3)

%Constantes de Adsorption a temperatura T
K1=(2.10E-2)*exp(-33.4E3/(R*(T)));
K4=(2.40E-2)*exp(-31.5E3/(R*(T)));

%Constantes de equilibrio
Keqr=10^((2073/T)-2.029);
Keq=1/Keqr;

% Velocidades de reaccion son funcion de las variables x(fracciones molares) y el tiempo
r1=@(t)k1(t)*(((x(1)*x(2))-(x(3)*x(4)/Keq))/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^2);
r2=@(t)k2(t)*((x(3)*(x(2)^0.2))/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.2);
r3=@(t)k3(t)*((x(3)*(x(2)^0.2))/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.2);
r4=@(t)k4(t)*((x(3)*(x(2)^0.2))/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.2);
r5=@(t)k5(t)*((x(3)*(x(2)^0.2))/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.2);
r6=@(t)k6(t)*((x(3)*(x(2)^0.2))/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.2);
r7=@(t)k7(t)*((x(3)*(x(2)^0.2))/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.2);
r8=@(t)k8(t)*((x(3)*(x(2)^0.2))/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.2);
r9=@(t)k9(t)*((x(3)*(x(2)^0.2))/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.2);
r10=@(t)k10(t)*((x(6)^1.3)/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.3);
r11=@(t)k11(t)*((x(7)^1.8)/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.8);
r12=@(t)k12(t)*((x(9)^1.6)/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.6);
r13=@(t)k13(t)*((x(10)^1.3)/(1+(K4*x(4)*pt)+(K1*x(1)*pt)))*(pt^1.3);

%Condiciones de alimentacion y calculo de la fraccion en masa de co2
nco2=1;
nh2=3;
nt=nco2+nh2;
mt=(nco2/nt*(M1)+nh2/nt*(M2));
wco2=(nco2*(M1)/mt);

%Faktor
% F=@(x) ((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2

%Funcion de las ecuaciones diferenciales en funcion de 'x' y 'r'
fg=@(t,x) ([ ((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*-r1(t);...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(-r1(t)-3*r2(t)-2*r3(t)-2*r4(t)-7/3*r5(t)-2*r6(t)-2*r7(t)-2*r8(t)-6/4*r9(t));...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(r1(t)-r2(t)-r3(t)-r4(t)-r5(t)-r6(t)-r7(t)-r8(t)-r9(t));...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(r1(t)+r2(t)+r3(t)+r4(t)+r5(t)+r6(t)+r7(t)+r8(t)+(2/4*r9(t)));...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*r2(t);...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(1/2*r3(t)-r10(t));...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(1/3*r4(t)-r11(t));...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(1/3*r5(t));...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(1/4*r6(t)-r12(t));...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(1/5*r7(t)-r13(t));...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(1/10*r8(t)+1/5*r10(t)+1/3.3*r11(t)+1/2.5*r12(t)+1/2*r13(t));...
((M1*x(1)+M2*x(2)+M3*x(3)+M4*x(4)+M5*x(5)+M6*x(6)+M7*x(7)+M8*x(8)+M9*x(9)+M10*x(10)+M11*x(11)+M12*x(12)))*wco2*(1/4*r9(t))]);

%Solucion
[t,x]=ode113(fg,[0,tf],x0);
% plot(t,x(:,1) )
figure(gcf)
axis([0 1 0 1 ])
hold on
for i=1:size(x,2)
plot(t,x(:,i),'color', rand(1,3) ,'linewidth',2 )
a{i}=['X=',num2str(i)];
legend(a)
axis([0 1 0 1])
grid on
xlabel('t')
ylabel('x , y');
title('Ben-Gurion Model')
pause(1)
end

__________________________________________________________________________________________

De hecho si es que cambio el vector de valores iniciales a:
https://www.lawebdelprogramador.com/img/img.png?11.4
x0=[0.05,0.55,0.2,0.2,0,0,0,0,0,0,0,0]. aqui les envio una foto con el nombre: BgV1(0.05-0.55-0.20-0.20)

el sistema se comporta mejor--> pero no se pq no funciona la generacion de x5 hasta x12 a partir de la condicion iniciales reales de:

x0=[0.25,0.75,0,0,0,0,0,0,0,0,0,0]

ya que como les dije el sistema si genera x(3). Segun el modelo matematico me imagino el siguiente comportamiento de la variable x3 (llegando a un maximo) les envio aqui una foto como es que me imagino el comportamiento de las variables segun el sistema de ecuaciones que implemente. El nombre de la segunda foto es:

BgV1-Vstlg

Gracias de antemano

• BgV10.05-0.55-0.20-0.20.zip(94,5 KB)

Hay preguntas que necesitan tiempo para poder entender porque hay vacíos en las repreguntas planteadas y lo recomendable es una asesoría personalizada, donde el usuario que hace la pregunta explique en detalle lo que desea y desde ese punto de vista se puede hacer un código de acuerdo a su necesidad o en su defecto corregir el código ya avanzado por el usuario que hace la pregunta.

Sugerencia:
Trate de seguir haciendo las preguntas en la pregunta original que hizo, para que cualquier usuario de matlab que tenga ese sentir de ayudarle, pueda entenderder su pregunta y pueda apoyarlo.

https://www.lawebdelprogramador.com/foros/Matlab/1683035-ode113-y-cinetica-quimica-Undefined-function-or-variable-x.html

Saludos
JOSE JEREMIAS CABALLERO
Asesor de Proyectos con Matlab
Servicios de programación matlab

http://matlabcaballero.blogspot.com
https://www.facebook.com/matlabcaballero

Muchas gracias

Estimado Sr. Caballero

el codigo estaba bien y asi tiene que funcionar, fue un error de concepto, el que estaba cometiendo ay ay ay... Le agradezco por su tiempo y seguramente hasta en breve saludos :)