ode45 not working (giving non real values)

Why am I getting non real values for M.

iFunction file

function dMdx=s(x,M,gamma_b,a1Vec,b1Vec,a11Vec,b11Vec,xVec,xMin)
a1=interp1(xVec,a1Vec,x);
b1=interp1(xVec,b1Vec,x);
a11=interp1(xVec,a11Vec,x);
b11=interp1(xVec,b11Vec,x);
if xVec==xMin
    dMdx=(1/4)*(-gamma_b*a1+sqrt((gamma_b*a1)^2-4*((gamma_b-1)*a1^2-(1+gamma_b)*a11+(((1+gamma_b)^2)/2)*b11)));
else
dMdx= M*((1+((gamma_b-1)/2)*M^2)/(1-M^2))*(-a1+((1+(gamma_b*M^2))/2)*b1);
end
end

iScript

xVec=0:0.01:0.4;  %axial combuster length (can be varied according to user input and units)
xi=xVec(1);  %station 3 (combuster starting point)
x4=xVec(end);%station 4 (combuster end point)
gamma_b=1.36;%Average value of ratio of specific heat capacities during burning
gamma_c=1.4;%Average value of ratio of specific heat capacities during compression
X=((xVec-xi)/(x4-xi));
a3=0.0038;%station 3 cross secion area in m2.
A=a3*(1+X.^2);%Area Profile
a=gradient(A,xVec);%dA/dx in equation (6-90)
a1Vec=a./A; %(dA/dx)/A in equation (6-90)
d2a=gradient(a,xVec);
a11Vec=d2a./A;
tau_b=1.5; %Tt4/Tt2 in eq. (6-91) overall total temerature rise in burner
theta=2;%emperical constant in eq. (6-91)
Tt2=2200;%Total temperature at station 2 depends on inlet and compression conditions (USER SPECIFIED)
Ttx=Tt2*(1+(tau_b-1)*((theta*X)./(1+(theta-1)*X)));%Temperature profile eq. (6-91)
b=gradient(Ttx,xVec);%dTt/dx in equation (6-90)
b1Vec=b./Ttx;%(dTt/dx)/Tt in equation (6-90)
d2b=gradient(b,xVec);
b11Vec=d2b./Ttx;
fun = @(xVec)((a./A)-(((1+gamma_b)/2)*(b./Ttx))); % Function of x only
[fMin, idx] = min(abs(fun(xVec)));
xMin = xVec(idx);
initial_M=1;
[x,M]=ode45(@(x,M) s(x,M,gamma_b,a1Vec,b1Vec,a11Vec,b11Vec,xVec,xMin),[xMin xi],initial_M);%Solve ODE eq. (6-90)

Answer 1

[Walter Roberson]

Your test

if xVec == xMin

is testing a vector of values against a single value. Your xVec are distinct so at most once of those will compare equal, so you will get a logical vector that has at most one true. When you have an "if" that is being asked to check a vector, the result is considered true only if all of the values in the vector are true. That can never be the case for you so the "else" is always taken.

You probably want to test x == xMin instead of xVec == xMin. But if you do that you should probably keep in mind that it is not easy to establish conditions under which floating point numbers that are intended to be the same value will be bit-for-bit identical and so end up comparing true. Your code as it stands does use the right conditions so if you do not change it then if you compare x == xMin then it will compare true the first time through, but it is probably not a good idea to count on that.