Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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
share|improve this question
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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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