# Solving a system of nonlinear coupled ODEs

I tried using NDSolve to solve the following system of equations:

t x'[t] == -x[t] + y[t],
t y'[t] == -5 t^2/x[t]^2 + x[t] - y[t],
x[1] == 4, x[100] == 1


It's weird. The system tells me there is infinity at the boundary $t=1$, so I change the boundary from $t=1$ to $t=2$, and I get the same message again.

If I eliminate the nonlinear term -5 t^2/x[t]^2, the function can be solved analytically. So I do not know whether these equations are well posed or not with such boundary conditions just for x[t]

Anyone have a suggestion?

• Please post the precise command you entered and quote the error message you got. Commented Sep 28, 2013 at 23:28
• putting a side for a minute the 1/0 error, you need to have an initial/boundary condition for $y(t)$. Even though these are coupled, each derivative in the equations produces one constant of integration on its own. Are you sure you copied this problem correctly? Commented Sep 29, 2013 at 2:50
• @Nasser, yes, I do need to specify a boundary condition for y(t) to satisfy my boundary conditions for x(t)
– 3c.
Commented Sep 29, 2013 at 17:31

As noted above, you need an initial condition for y. You can then solve the odes given x[1] == 4 and y[1] == y1 where y1 is some real number. By solving the equations for different values of y1, you can see if there is a value such that x[100] == 1. In Version 9, it turns out the ParametricNDSolve will help you do this. Make y1 a parameter:

{xsoln, ysoln} = {x, y} /. ParametricNDSolve[{
t x'[t] == -x[t] + y[t],
t y'[t] == -5 t^2/x[t]^2 + x[t] - y[t],
x[1] == 4,
y[1] == y1
}, {x, y}, {t, 1, 100}, {y1}]


Now locate a region where y1 is likely to give the desired result:

Plot[xsoln[y1][100], {y1, 43.35, 43.5}, GridLines -> {{}, {1}}]


Next use FindRoot to get a precise value and call it y1star:

y1star = y1 /. FindRoot[xsoln[y1][100] == 1, {y1, 43.4}]


Finally, confirm that it produces the desired result

Plot[{xsoln[y1star][t], ysoln[y1star][t]}, {t, 1, 100},
PlotRange -> {-10, All}, GridLines -> {{}, {1}}]


• thanks. I do not know whether my equations are well posed or not. you provide me a method to check whether my boundary conditions just for x[t] can work or not
– 3c.
Commented Sep 29, 2013 at 17:29

Although the solution by Nasser is entirely satisfactory, this alternative approach may be of interest. Designate the two equations in the question as eq1 and eq2, respectively. Then, since y has no boundary conditions, eliminate it from the system.

Flatten@Solve[Subtract @@ eq1 + Subtract @@ eq2 == 0, y'[t]]
eq3 = Simplify[D[eq1, t] /. %]
(* {y'[t] -> (-5 t - x[t]^2 x'[t])/x[t]^2} *)
(* (5 t)/x[t]^2 + 3 x'[t] + t x''[t] == 0 *)


Then, solve the resulting equation for x, and from that result determine y as well (from eq1).

{sx, sxp} = NDSolveValue[{eq3, x[1] == 4, x[100] == 1}, {x, x'}, {t, 1, 100},
Method -> {"Shooting", "StartingInitialConditions" -> x'[1] == 39.5}];
Plot[{sx[t], sx[t] + t sxp[t]}, {t, 1, 100}, PlotRange -> {-20, All}]


As pointed out by Nasser, you need to have an initial/boundary condition for y(t). Here I assumed a random condition on y(t)

sol = NDSolve[{t x'[t] == -x[t] + y[t], t y'[t] == -5 t^2/x[t]^2 + x[t] - y[t],
x[1] == 4, y[1] == 1}, {x, y}, {t, 1, 2}]

Plot[{x[t] /. sol, y[t] /. sol}, {t, 1, 2}, PlotRange -> All,
PlotStyle -> {Blue, Red}, PlotLegends -> Placed[{"x[t]", "y[t]"}, Above]]