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[0] == To,
   q[0] == 0.,
   nCO[0] == nCOinit,
   nO[0] == 0.,
   nO2[0] == nO2init,
   nCO2[0] == 0.,
   m[0] == nCOinit + nO2init
  {T, q, Tenv, nCO, nO, nO2, nCO2, m},
  {t, 0.00, tfinal},
  MaxStepSize -> 0.001,
  MaxSteps -> Infinity
  • 2
    $\begingroup$ 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? $\endgroup$
    – rcollyer
    Mar 11, 2012 at 2:53
  • 2
    $\begingroup$ 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. $\endgroup$
    – David
    Mar 11, 2012 at 3:15

1 Answer 1


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.