MATLAB: Does MATLAB return a solution with all “NaN”s and throw the warning message “Matrix is singular to working precision.” when using ODE15s

diffusionMATLABode'srcondsingularthermodynamics

I am trying to integrate a differential equation with a boundary condition that varies with time. For example, the call to "ode15s" below outputs a vector "y" full of "NaN"s:
>>odefun=@(T,Y)Ydot_TNES1(T,Y,NSTEP,TEMPK0,CR,MESHSP,DISLOC,DISPV,DISPI);
>>[T,Y]=ode15s(odefun,tspan,Y0,options);
What is the cause of this behavior?

Best Answer

In this case, you have a condition defined in the MATLAB file "Ydot_TNES1" that is changing with time. It is given by the expressions below. Note that "EFV" and "EFI" are constants, while "BOLTZ" is Boltzmann's constant:
>>CVDOT(1,1)=(EFV/BOLTZ*exp(-EFV/(BOLTZ*T)))/T^2;
>>CIDOT(1,1)=(EFI/BOLTZ*exp(-EFI/(BOLTZ*T)))/T^2;
After executing the code, notice that MATLAB throws the warning "Matrix is singular to working precision." This indicates that the problem is either poorly scaled, or singular. 
You can resolve the issue by noting that the variable "T" in the equations above is interpreted as "time" by ode15s. This can be verified by examining the function declaration syntax, as all of the ODE Suite functions expect function handles of the form "@(t,y)f(t,y)". At time "T = 0", both of the equations above are singular, so integration forward in time will not yield a finite solution. 
This issue is resolved by noting that this boundary condition is erroneous, and is due to a confusion of the variables in the problem. "T" is supposed to correspond to temperature. Upon changing "T" to "TEMP", the code runs without throwing any error messages, and the resulting solution is well scaled.
Note that, when debugging issues related to thermodynamics, it is important to remember that "Temperature" and "Time" are easily confused. One way to spot this issue quickly is to look for equations that have the expression below somewhere:
>>exp(-E/(BOLTZ*T));
This equation corresponds to the famous "Boltzmann Distribution", and "T" should be interpreted as temperature here.