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 want to solve the system of non-linear differential equations given below numerically.

$$y''(t)+500y'(t)+100y(t)=-33\cos(500t)-66\cos(1000t)$$ $$300x'(t)=1000y(t)+500y'(t)-35\tanh(50x'(t))$$

Notice that there is a derivative of $x$ inside $\tanh$. I tried the code below, but many errors appeared.

sol1 = 
  NDSolve[{y''[t] + 500 y'[t] + 100 y[t] == -33 Cos[500 t] - 66 Cos[1000 t],
    y[0] == 0, y'[0] == 0}, y, {t, 0, 30}];
sol2 = 
  NDSolve[{300 x'[t] == 
    1000 First[Evaluate[y[t] /. sol1]] + 500 First[Evaluate[y'[t] /. sol1]] - 
    35 Tanh[50 x'[t]], x[0] == 0}, x, {t, 0, 30}];
Plot[{Evaluate[x[t] /. sol1], Evaluate[y[t] /. sol2]}, {t, 0, 30}, PlotRange -> All]

How can I solve my system of differential equations using Mathematica?


I found the answer! thanks to Nasser. The solution for $x(t)$ is like below. I just need to adjust the axis

enter image description here

share|improve this question
up vote 2 down vote accepted

Is the solution supposed to be something like this? Hopefully I did not break too many things. You should know if the solution looks right or not.

Clear[y, x, t];
eqs = y''[t] + 500 y'[t] + 100 y[t] == -33 Cos[500 t] - 66 Cos[1000 t];
ic = {y[0] == 0, y'[0] == 0};
sol1 = First@NDSolve[{eqs, ic},{y[t], y'[t]}, {t, 0, 30},Method -> "StiffnessSwitching"];
eqs = 300 x'[t] == 1000 Evaluate[y[t] /. sol1] + 500 Evaluate[y'[t] /. sol1] - 
   35 Tanh[50 x'[t]]

Mathematica graphics

sol2 = NDSolve[{eqs, x[0] == 0}, x[t], {t, 0, 30},MaxSteps -> Infinity]

Mathematica graphics

Grid[{{Plot[Evaluate[y[t] /. sol1], {t, 0, 0.5}, PlotRange -> All]},
  {Plot[Evaluate[x[t] /. sol2], {t, 0, 30}, PlotRange -> All]}}]

Mathematica graphics

share|improve this answer
thank you! i edited my question to include the answer. – Dadan Ari Wibowo Nov 11 '13 at 6:50
But if it is okay, I would like to know the reason for using StiffnessSwitching and MaxSteps->Infinity... – Dadan Ari Wibowo Nov 11 '13 at 6:56
"StiffnessSwitching", well, more information on this is here…;. Physics wise, the natural frequency of your eqs is 10 rad/sec, while the forcing frequencies (RHS) are much higher (you have one at 1000 rad/sec). So it tells one that system will oscillate slowly in the early transient stage, then in steady state, when the load takes over, system will oscillate at much higher frequency (at the driving frequency). Hence you have a system with low and high vibrations in it. Hence this method to try.... – Nasser Nov 11 '13 at 7:10
..."MaxSteps->Infinity" this one is easy. Mathematica gave an error saying Maximum number of 10000 steps reached at the point t == 1.5727575203, so you just increase the number of steps and hope for the best. – Nasser Nov 11 '13 at 7:11

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.