Hola, Gracias por la respuesta, en el transcurso del sía continue analizando ara encontrar el error y lo encontré, ahora ya no me dice que la matriz es singular pero aparece "Busy" y se queda asi, lo dejé 1 hora corriendo y continuaba así, colocaré las ecuaciones tal como las escribé en matlab.
Los "W" representan los ratios de comsumo o producción de las especies (derivadas de la concentración en función del tiempo), y todas las otras variables que aparecen son funciones de la temperatura, la función dy(15) representa la derivada de la temperatura en función del tiempo.
Por ejemplo
kf1 = A(1)*(y(15))^n(1)*exp(-E(1)/(R1*y(15))), las otras kf son funciones del mismo tipo con valores diferentes de A,n y E.
kb son determinadas a partir de las kf utilizando las energías de gigbs de la reacción. Las energías de Gigbs son polinomios de la temperatura.
Las capacidades caloríficas (Cp) y las entalpías (H) son polinomios de la temperatura.
WO2 = -kf1*y(1)*y(4)+kb1*y(5)*y(3)+kf6*y(3)*y(3)*(y(10)+2.4*y(2)+15.4*y(6)+1.75*y(8)+3.6*y(9))-kb6*y(1)*(y(10)+2.4*y(2)+15.4*y(6)+1.75*y(8)+3.6*y(9))-kf11*y(1)*y(8)+kb11*y(9)*y(3)-kf13*y(1)*y(11)+kb13*y(12)*y(3)+kf17*y(14)*y(3)-kb17*y(12)*y(1)-kf18*y(1)*y(4)*(y(10)+0.75*y(8)+1.5*y(9))+kb18*y(7)*(y(10)+0.75*y(8)+1.5*y(9))+kf21*y(13)*y(3)-kb21*y(10)*y(1);
WH2 = -kf2*y(2)*y(3)+kb2*y(5)*y(4)-kf3*y(2)*y(5)+kb3*y(6)*y(4)+kf5*y(4)*y(4)*y(10)-kb5*y(2)*y(10);
WO = kf1*y(1)*y(4)-kb1*y(5)*y(3)-kf2*y(2)*y(3)+kb2*y(5)*y(4)+kf4*y(5)*y(5)-kb4*y(6)*y(3)-kf6*y(3)*y(3)*(y(10)+2.4*y(2)+15.4*y(6)+1.75*y(8)+3.6*y(9))+kb6*y(1)*(y(10)+2.4*y(2)+15.4*y(6)+1.75*y(8)+3.6*y(9))-kf10*y(8)*y(3)*(y(10)+2*y(2)+6*y(1)+6*y(6)+1.5*y(8)+3.5*y(9))+kb10*y(9)*(y(10)+2*y(2)+6*y(1)+6*y(6)+1.5*y(8)+3.5*y(9))+kf11*y(1)*y(8)-kb11*y(9)*y(3)+kf12*y(12)*y(11)-kb12*y(3)*y(10)+kf13*y(1)*y(11)-kb13*y(12)*y(3)-kf17*y(14)*y(3)+kb17*y(11)*y(1)+kf19*y(13)*(y(10)+2*y(2)+6*y(6)+1.5*y(8)+2*y(9))-kb19*y(10)*y(3)*(y(10)+2*y(2)+6*y(6)+1.5*y(8)+2*y(9))-kf20*y(13)*y(3)+kb20*y(11)*y(11)-kf21*y(13)*y(3)+kb21*y(10)*y(1);
WH = -kf1*y(1)*y(4)+kb1*y(5)*y(3)+kf2*y(2)*y(3)-kb2*y(5)*y(4)+kf3*y(2)*y(5)-kb3*y(6)*y(4)-kf5*y(4)*y(4)*y(10)+kb5*y(2)*y(10)-kf7*y(4)*y(5)*(y(10)+0.73*y(2)+3.65*y(6))+kb7*y(6)*(y(10)+0.73*y(2)+3.65*y(6))+kf8*y(8)*y(5)-kb8*y(9)*y(4)+kf14*y(11)*y(5)-kb14*y(12)*y(4)-kf16*y(14)*y(4)+kb16*y(12)*y(5)-kf18*y(1)*y(4)*(y(10)+0.75*y(8)+1.5*y(9))+kb18*y(7)*(y(10)+0.75*y(8)+1.5*y(9));
WOH = kf1*y(1)*y(4)-kb1*y(5)*y(3)+kf2*y(2)*y(3)-kb2*y(5)*y(4)-kf3*y(2)*y(5)+kb3*y(6)*y(4)-kf4*y(5)*y(5)+kb4*y(6)*y(3)-kf7*y(4)*y(5)*(y(10)+0.73*y(2)+3.65*y(6))+kb7*y(6)*(y(10)+0.73*y(2)+3.65*y(6))-kf8*y(8)*y(5)+kb8*y(9)*y(4)+kf9*y(8)*y(7)-kb9*y(9)*y(5)-kf14*y(11)*y(5)+kb14*y(12)*y(4)+kf15*y(12)*y(7)-kb15*y(14)*y(5)+kf16*y(14)*y(4)-kb16*y(12)*y(5);
WH2O = kf3*y(2)*y(5)-kb3*y(6)*y(4)+kf4*y(5)*y(5)-kb4*y(6)*y(4)+kf7*y(4)*y(5)*(y(10)+0.73*y(2)+3.65*y(6))-kb7*y(6)*(y(10)+0.73*y(2)+3.65*y(6));
WHO2 = -kf9*y(8)*y(7)+kb9*y(9)*y(5)+kf18*y(1)*y(4)*(y(10)+0.75*y(8)+1.5*y(9))-kb18*y(7)*(y(10)+0.75*y(8)+1.5*y(9));
WCO = -kf8*y(8)*y(5)+kb8*y(9)*y(4)-kf9*y(8)*y(7)+kb9*y(9)*y(5)-kf10*y(8)*y(3)*(y(10)+2*y(2)+6*y(1)+6*y(6)+1.5*y(8)+3.5*y(9))+kb10*y(9)*(y(10)+2*y(2)+6*y(1)+6*y(6)+1.5*y(8)+3.5*y(9))-kf11*y(8)*y(1)+kb11*y(9)*y(3);
WCO2 = kf8*y(8)*y(5)-kb8*y(9)*y(4)+kf9*y(8)*y(7)-kb9*y(9)*y(5)+kf10*y(8)*y(3)*(y(10)+2*y(2)+6*y(1)+6*y(6)+1.5*y(8)+3.5*y(9))-kb10*y(9)*(y(10)+2*y(2)+6*y(1)+6*y(6)+1.5*y(8)+3.5*y(9))+kf11*y(8)*y(1)-kb11*y(9)*y(3);
WN2 = kf12*y(12)*y(11)-kb12*y(3)*y(10)+kf19*y(13)*(y(10)+2*y(2)+6*y(6)+1.5*y(8)+2*y(9))-kb19*y(3)*y(10)*(y(10)+2*y(2)+6*y(6)+1.5*y(8)+2*y(9))+kf21*y(3)*y(13)-kb21*y(1)*y(10);
WN = -kf12*y(12)*y(11)+kb12*y(3)*y(10)-kf13*y(11)*y(1)+kb13*y(11)*y(3)-kf14*y(11)*y(5)+kb14*y(12)*y(4);
WNO = -kf12*y(12)*y(11)+kb12*y(3)*y(10)+kf13*y(11)*y(1)-kb13*y(12)*y(3)+kf14*y(11)*y(5)-kb14*y(12)*y(4)-kf15*y(12)*y(7)+kb15*y(14)*y(5)+kf16*y(14)*y(4)-kb16*y(12)*y(5)+kf17*y(14)*y(3)-kb17*y(12)*y(1)+kf20*y(13)*y(3)-kb20*y(12)*y(12);
WN2O = -kf19*y(13)*(y(10)+2*y(2)+6*y(6)+1.5*y(8)+2*y(9))+kb19*y(10)*y(3)*(y(10)+2*y(2)+6*y(6)+1.5*y(8)+2*y(9))-kf20*y(13)*y(3)+kb20*y(12)*y(12)-kf21*y(13)*y(3)+kb21*y(10)*y(1);
WNO2 = kf15*y(12)*y(7)-kb15*y(14)*y(5)-kf16*y(14)*y(4)+kb16*y(12)*y(5)-kf17*y(14)*y(3)+kb17*y(12)*y(1);
%Differenctial Equations
dy=zeros(15,1);
dy(1) = WO2;
dy(2) = WH2;
dy(3) = WO;
dy(4) = WH;
dy(5) = WOH;
dy(6) = WH2O;
dy(7) = WHO2;
dy(8) = WCO;
dy(9) = WCO2;
dy(10)= WN2;
dy(11)= WN;
dy(12)= WNO;
dy(13)= WN2O;
dy(14)= WNO2;
dy(15)= (R*y(15)*(WO2+WH2+WO+WH+WOH+WH2O+WHO2+WCO+WCO2+WN2+WN+WNO+WN2O+WNO2)-(HO2*WO2+HH2*WH2+HO*WO+HH*WH+HOH*WOH+HH2O*WH2O+HHO2*WHO2+HCO*WCO+HCO2*WCO2+HN2*WN2+HN*WN+HNO*WNO+HN2O*WN2O+HNO2*WNO2))/(y(1)*(Cp_O2-R)+y(2)*(Cp_H2-R)+y(3)*(Cp_O-R)+y(4)*(Cp_H-R)+y(5)*(Cp_OH-R)+y(6)*(Cp_H2O-R)+y(7)*(Cp_HO2-R)+y(8)*(Cp_CO-R)+y(9)*(Cp_CO2-R)+y(10)*(Cp_N2-R)+y(11)*(Cp_N-R)+y(12)*(Cp_NO-R)+y(13)*(Cp_N2O-R)+y(14)*(Cp_NO2-R));
Bueno como dije ahora el programa se toma demasiado tiempo y no me da una respuesta, no sé si en algún momento lo haría suponiendo que lo deje correr por mucho más. Como vez la forma en que escribí las ecuaciones no es muy sofisticada pero no sé si eso justifica que el programa tarde tanto.
Gracias por la respuesta.