# 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 == 4, x == 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. – Szabolcs Sep 28 '13 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? – Nasser Sep 29 '13 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. Sep 29 '13 at 17:31

As noted above, you need an initial condition for y. You can then solve the odes given x == 4 and y == 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 == 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 == 4,
y == y1
}, {x, y}, {t, 1, 100}, {y1}]


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

Plot[xsoln[y1], {y1, 43.35, 43.5}, GridLines -> {{}, {1}}] Next use FindRoot to get a precise value and call it y1star:

y1star = y1 /. FindRoot[xsoln[y1] == 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. Sep 29 '13 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 == 4, x == 1}, {x, x'}, {t, 1, 100},
Method -> {"Shooting", "StartingInitialConditions" -> x' == 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 == 4, y == 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]] 