Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm having trouble using NDSolve to solve equations that were obtained from solving cubic equations. In my problem, I first obtained 6 functions (ca1, cb1,cc1, ca2,cb2,cc2) that are in terms of a1[t], b1[t],c1[t], a2[t],b2[t],c2[t]. So it will be like this:







So what I did is :

Solve[ca1 == k1*b1[t]*c1[t]^3*(1- 3 ca1 - 2 cb1 - cc1)^3 && 
  cb1 == k2*b1[t]*c1[t]^4*(1- 3 ca1 - 2 cb1 - cc1)^2 && 
  cc1 == k3*c1[t]^2*(1- 3 ca1 - 2 cb1 - cc1), {ca1, cb1, cc1}]

Solve[ca2 == k1*b2[t]*c2[t]^3*(1- 3 ca2 - 2 cb2 - cc2)^3 && 
      cb2 == k2*b2[t]*c2[t]^4*(1- 3 ca2 - 2 cb2 - cc2)^2 && 
      cc2 == k3*c2[t]^2*(1- 3 ca2 - 2 cb2 - cc2), {ca2, cb2, cc2}]

And then I obtained a very long solution (that contains complex solution) for ca1, cb1,cc1, ca2,cb2,cc2 . And I tried to used the 'uncomplex solution' that does not contain any complex number, 'i'. And then I continued to construct the final set of the equations:

D1 = 1.*10^-3;
D2 = 1.*10^-3;
D3 = 1.*10^-3;

K1 = .05;
K2 = 0.2;
K3 = .01;

e1 = {a1'[t] == D1*(ca1 - ca2),
   a2'[t] == -D1*(ca1 - ca2),
   b1'[t] == D2*(cb1 - cb2),
   b2'[t] == -D2*(cb1 - cb2),
   c1'[t] == D3*(cc1 - cc2),
   c2'[t] == -D3*(cc1 - cc2),
   a1[0] == 0.02, a2[0] == 10^-10, b1[0] == 4, b2[0] == 0.001, 
   c1[0] == 0.0001, c2[0] == 10^-10};

And just a quick note that ca1, ca2, cb1, cb2, cc1, and cc2 are functions of a1[t], a2[t], b1[t], b2[t], c1[t], c2[t].

Then applied NDSolve:

sol1 = NDSolve[e1, {a1, a2, b1, b2, c1, c2}, {t, 0, 3000}];

I tried to solve this problem on my Mathematica 8.0, but I kept getting a message saying " the max step reach at 1000". I looked up an older post on this forum saying we could used this code.

sol1 = NDSolve[e1, {a1, a2, b1, b2, c1, c2}, {t, 0, 3000}, 
   Method -> {"FixedStep", Method -> "ExplicitEuler"}, 
   StartingStepSize -> 10000];

But my solution looks funky...Could someone help telling me what I did wrong here? Or how I could solve this problem in better way? Thanks for you time reading this. I hope my explanation is clear...Otherwise, please let me know if I forgot to provide any information.

share|improve this question
I think you want to set MaxSteps->Infinity instead of StartingStepSize and remove the other method options. –  user21 Feb 21 '13 at 18:39
@ ruebenko, I tried 'MaxSteps->Infinity'. But it seems to take a really long time to calculate. ( It's been running for 15 min now but still calculating). Thanks. –  DumbleKo Feb 21 '13 at 19:03
you can add a Dynamic["time: " <> ToString[CForm[currentTime]] <> " step: " <> ToString[CForm[steps]]] steps = 0; before the call to NDSolve and then add the following option into the call of NDSolve: EvaluationMonitor :> (currentTime = t; steps++). Does it get stuck at a point or is it 'just' slow? –  user21 Feb 21 '13 at 23:43
@ ruebenko, your method in your first reply work on all my cases that does not involved cubic equation. However, I still get error message from using your second reply. Thanks. –  DumbleKo Feb 24 '13 at 5:12
@DunbleKo, when I copy your code and eval, it does not give the error you claim; can you please fix this, otherwise I have no chance to help. –  user21 Feb 24 '13 at 9:11
show 1 more comment

Your Answer


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

Browse other questions tagged or ask your own question.