ode return NaN :(

clc;clear;
Gamma = 5/3;
Lambda = (Gamma+1)/2;
syms y(t) x(t) 
eqn1 = diff(y,2)*(y-x) == Lambda/2*diff(y)*(diff(x)-diff(y)/Lambda);
a = sqrt(Gamma*(Lambda-1))/Lambda;
dts = (y-x)/(a*diff(y));
Pp = 1/0.08*(t-dts);
pip = Pp/(1/Lambda);

eqn2 = diff(x,2)*y == (pip-diff(y)^2/Lambda);
eqn = [eqn1 eqn2];
[V,S] = odeToVectorField(eqn);
M = matlabFunction(V,'vars',{'t','Y'});
interval = [0 2];
yInit = [0 0 0 0];
ySol = ode45(M,interval,yInit);
tValues =  linspace(0,2,10);
xValues = deval(ySol,tValues,1);
zValues = deval(ySol,tValues,3);

xValues =

0 NaN NaN NaN NaN NaN NaN NaN NaN NaN

In fact, I want to calculate this equation.

[Torsten]

You divide by Zdot which is 0 for t=0.

Answer 1

[Sam Chak]

Hi @ZH CC

The odeToVectorField(eqn) returns this

V =

Y[2]

-(9*Y[4]^3 - 160*5^(1/2)*Y[1] + 160*5^(1/2)*Y[3] - 200*t*Y[4])/(12*Y[3]*Y[4])

Y[4]

-((4*Y[2] - 3*Y[4])*Y[4])/(6*(Y[1] - Y[3]))

S =

x

Dx

y

Dy

and you can clearly see the divisions by and . Since the initial values are , and , naturally, the ode45 solver returns NaN.

Hope you are satisfied with this Answer.

[ZH CC]

Thanks a lot, this is really good tips.

Answer 2

[Steven Lord]

What is at time t = 0? By your formula 14 it is times other stuff. Let's ignore that other stuff and look at the value of that term. Well, all of Z(0), X(0), and (0) are 0 by your formula 17. So that term is which MATLAB correctly computes as NaN.

My guess is that at least one of Z(0), X(0), and (0) must not be 0. In particular, because Z-X appears in the denominator I think one or both of those must be non-zero (and not equal to each other) for your equations to make sense.

[ZH CC]

Thank you very much for your answer.