ode45 returning Nan values

Iam trying to solve a system of differential equations using the ode45 command. However I am only getting Nan values as the output. Here's the code I have so far-

clc
global rho Cp MW dH

tspan = 60*[1 300000];
y0 = [1 333.15];
[t, y] = ode45(@myODE, tspan, y0)
plot(t/60, y(:,2)-273.15)

function dydt = myODE(t, y)
    rho = 906;
    Cp = 0.4365;
    MW = 104.15;
    dH = -17800;
    
    y2 = 1-y(1);
    
    A0 = 1.964*(10^5)*exp(-10040/y(2));
    A1 = 2.57-5.05*y(2)*(10^(-3));
    A2 = 9.56-1.76*y(2)*(10^(-2));
    A3 = -3.03+7.85*y(2)*(10^(-3));
    A = A0*exp(A1*(y2) + A2*(y2^2) + A3*(y2^3));
    
    dydt = zeros(2,1);
    dydt(1) = -A*((rho*y(1)/MW)^(1.5))*y(1);
    dydt(2) = (dH/(Cp*MW))*(-A*((rho*y(1)/MW)^(1.5))*y(1));
end

This is the output I get-

y =

1.0e+02 *

0.0100 + 0.0000i 3.3315 + 0.0000i

NaN + NaNi NaN + NaNi

NaN + NaNi NaN + NaNi

NaN + NaNi NaN + NaNi

I am not sure where I am going wrong here. Is there a problem with the equations? Any help is appreciated!

Answer 1

[Alan Stevens]

Have you basically got the wrong sign in dydt(2)?

tspan = 60*[1 300000];
y0 = [1 333.15];
[t, y] = ode45(@myODE, tspan, y0);
plot(t/60, y(:,2)-273.15)

function dydt = myODE(~, y)
    rho = 906;
    Cp = 0.4365;
    MW = 104.15;
    dH = -17800;
    
    y2 = 1-y(1);
    
    A0 = 1.964*(10^5)*exp(-10040/y(2));
    A1 = 2.57-5.05*y(2)*(10^(-3));
    A2 = 9.56-1.76*y(2)*(10^(-2));
    A3 = -3.03+7.85*y(2)*(10^(-3));
    A = A0*exp(A1*(y2) + A2*(y2^2) + A3*(y2^3));
    
    dydt = zeros(2,1);
    dydt(1) = -A*((rho*y(1)/MW)^(1.5))*y(1);
    dydt(2) = (dH/(Cp*MW))*(A*((rho*y(1)/MW)^(1.5))*y(1));   %  Changed the sign from -A to +A
end

[Nishant Gupta]

No that equation is correct. The negative sign should be there

[Alan Stevens]

Ok. What about in dydt(1)?

Do you expect the temperatures to increase or decrease?

[Nishant Gupta]

I expect the temperature to increase. As dH is negative, heat will be evolved as the reaction proceeds. y(1) is the mole fraction and y(2) is the temp.