# How can I get results over the whole domain of my non-linear differential equations?

I want to solve the following system of non-linear differential equations:

t0 = 0.6;
tmax = 20;

σ = 0.04 t0;
a1 = -0.02;
a2 = -0.8;

s =
NDSolve[
{D[x[t], t] + σ x[t] + a1/t0 (y[t]/t)^(1/3) == 0,
D[y[t], t] - a2 x[t]/t == 0, x[t0] == y[t0] == 1},
{x, y}, {t, t0, tmax}]

Plot[x[t] /. s, {t, t0, tmax}]


The problem is that the NDSolve long before reaching max.

I guess the problem comes from the power 1/3.

How can I correct the situation?

• Please do not post images of your work. Please post your actual Mathematica code in the form of text that can be copied and pasted into a Mathematica notebook. Without such, it will be difficult to reproduce your problem and to experiment with possible solutions. – m_goldberg Feb 20 at 6:19

Surd is the solution, but it behaves delicately in this specific situation. Surd[x, n] returns the $$n$$-th real root of $$x$$, but it requires that $$n$$ be an integer. If NDSolve is fed any inexact coefficients, however, it will numericize everything -- including Surd's second argument it turns out.

Thus, use exact arithmetic to define the variables in the equation:

t0 = 6/10;
tmax = 20;
σ = 4/100 t0;
a1 = -2/100;
a2 = -8/10;


And then use Surd[y[t]/t, 3] in place of (y[t]/t)^(1/3):

s = NDSolve[{D[x[t], t] + σ x[t] + a1/t0 Surd[y[t]/t, 3] == 0,
D[y[t], t] - a2 x[t]/t == 0, x[t0] == 1, y[t0] == 1}, {x, y}, {t,
t0, tmax}]


This is also much faster to solve, likely because it doesn't encounter any abnormalities in the complex plane.

You can easily check that it matches up with the other solution (where the other solution is defined at all) by naming this solution something other than s and plotting both x[t] results on the same graph.