# NDSolve does not give me the expected solution

I have an assignment where I need to compute chemical reactions by solving ODEs. I use NDSolve to do this. However, for one problem it does not give me the predicted answer. My friend uses MATLAB's ode15s function to solve it and he gets the correct solution.

Is there any option in NDSolve that I should play with to match the correct solution?

Sorry if I don't give much information on the nature of the problem since that would be tedious. However, here is the code I am using:

sol[To_?NumericQ] := NDSolve[{
T'[t] == (-q[t]/V - (uCO[T[t]] nCO'[t] + uO[T[t]] nO'[t] + uO2[T[t]] nO2'[t] + uCO2[T[t]] nCO2'[t]) / ((CvCO[T[t]] nCO[t] + CvO[T[t]] nO[t] + CvO2[T[t]] nO2[t] + CvCO2[T[t]] nCO2[t]))),
q[t] == 4 Pi r Koil (T[t] - Tenv[t]),
Tenv[t] == To,
nCO'[t] == -nCO[t] nO[t] m[t] r1[T[t]] - nCO[t] nO2[t] r2[T[t]] + nCO2[t] m[t] r1b[T[t]] + nCO2[t] nO[t] r2b[T[t]],
nO'[t] == -nCO[t] nO[t] m[t] r1[T[t]] + nCO[t] nO2[t] r2[T[t]] - 2 nO[t] nO[t] m[t] r3[T[t]] + nCO2[t] m[t] r1b[T[t]] - nCO2[t] nO[t] r2b[T[t]] + 2 nO2[t] m[t] r3b[T[t]],
nO2'[t] == -nCO[t] nO2[t] r2[T[t]] + nO[t] nO[t] m[t] r3[T[t]] + nCO2[t] nO[t] r2b[T[t]] - nO2[t] m[t] r3b[T[t]],
nCO2'[t] == nCO[t] nO[t] m[t] r1[T[t]] + nCO[t] nO2[t] r2[T[t]] - nCO2[t] m[t] r1b[T[t]] - nCO2[t] nO[t] r2b[T[t]],
m'[t] == nCO'[t] + nO'[t] + nO2'[t] + nCO2'[t],
T == To,
q == 0.,
nCO == nCOinit,
nO == 0.,
nO2 == nO2init,
nCO2 == 0.,
m == nCOinit + nO2init
},
{T, q, Tenv, nCO, nO, nO2, nCO2, m},
{t, 0.00, tfinal},
MaxStepSize -> 0.001,
MaxSteps -> Infinity
]

• If it isn't giving you what you expect, then what is it giving you? In other words, can you post a plot of the correct solution and the incorrect one? – rcollyer Mar 11 '12 at 2:53
• Please provide a minimal example; if that's not possible (it usually is), at least make the long code you post executable on its own. – David Mar 11 '12 at 3:15

If your friend used MATLAB's ode15s solver, then your equations are most likely stiff, and the default NDSolve options are unlikely to give you a numerically stable result unless your step sizes are insanely small. You'll need a stiff solver instead.
For examples on how to proceed in Mathematica, look at tutorial/NDSolveStiffnessTest. Specifically, explore the methods "StiffnessSwitching", "StiffnessTest" and "NonStiffTest". Since you haven't described your system, there is nothing more I can say at this point, but these hints should give provide a good starting point for you.