# How to solve a pair of coupled 2nd order equations using shooting method

I have a pair of equations for x[t] and y[t]. Both are 2nd order and they are coupled. Now say that I know the following boundary conditions: the value of x[10], x[10000], y[10000] say. And I know how y'[10] is related to y[10]. Then how should I solve it?

I think the shooting method is required. Can anyone tell me how to implement shooting method in this case?

You don't need to write the whole routine for me. Just give me the basic idea. Everywhere I can see examples shown for using the shooting method for single differential equations.

Solution =
Map[First[NDSolve[x''[
t] + ((f'[t]/f[t]) + 2/t) x'[
t] + ((y[t]^2/f[t]^2) + (2/(f[t]))) x[t] ==
0, y''[t] + (2/t) y'[t] -
2 (x[t]^2/f[t])y[t] ==
0, x[t0 + .01] == (-3/2) t0 x'[r0 + .01], y[
t0 + .01] == 0, y[100000 + .01] ==
10, x[100000 + .01] == 0, {x[t],y[t]}, t,
Method -> {"Shooting",
"StartingInitialConditions" -> {x[
t0 + .01] == (-3/2) t0 x'[t0 + .01], x'[
t0 + .01] == #}}]]&,{1.5,3,3.5}]
f[t] := (t^2) - (M/t)
t0 = 2
M = 8


Thanks a lot.

• NDSolve[] can do shooting; barring that, you can implement shooting yourself with ParametricNDSolve[]. Without your equations, however, there's nothing much to say or do. Commented Aug 17, 2015 at 14:16
• We need some more details. A quick search of this site shows at least a few examples of the NDSolve being used with the shooting method. Could you try to rephrase your problem with more details, especially highlighting what doesn't work for you with the built-in method? Commented Aug 17, 2015 at 14:17
• Your code doesn't run as is. If you want a parametric solution, you might want to look into ParametricNDSolve as well. Commented Aug 17, 2015 at 14:46
• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Commented Aug 17, 2015 at 16:43

Basically, there are syntax errors. Fix those and it seems to work. In particular, read the tutorial/DefiningFunctions carefully. Also, all the equations of your differential equation problem should be enclosed in {}. Finally, you have a r0 which I interpreted to be a typo and read it as t0.

ClearAll[f];
f[t_] := (t^2) - (M/t);  (* <-- N.B. the underscore *)
t0 = 2;
M = 8;

Solution =
Map[(Print[#];
First[NDSolve[
(*     v--- Note the List {}    *)
{x''[t] + ((f'[t]/f[t]) + 2/t) x'[t] + ((y[t]^2/f[t]^2) + (2/(f[t]))) x[t] == 0,
y''[t] + (2/t) y'[t] - 2 (x[t]^2/f[t]) y[t] == 0,
x[t0 + .01] == (-3/2) t0 x'[t0 + .01],
y[t0 + .01] == 0,
y[100000 + .01] == 10,
x[100000 + .01] == 0},
{x, y}, {t, t0 + .01, 100000 + .01},
Method -> {"Shooting",
"StartingInitialConditions" -> {x[t0 + .01] == (-3/2) t0 x'[t0 + .01],
x'[t0 + .01] == #}}]]) &, {1.5, 3, 3.5}];
(*
1.5
3
3.5  --> Message below indicates that the third starting IC is not so good
*)


NDSolve::berr: The scaled boundary value residual error of 44.1719386003899 indicates that the boundary values are not satisfied to specified tolerances. Returning the best solution found. >>

Hold@Plot[{x[t], y[t]} // Evaluate, {t, t0 + .01, 10 + .01}] /.
Solution // ReleaseHold // GraphicsRow


The plot over the whole domain reveals only that a steady state is reached fairly soon.

Hold@Plot[{x[t], y[t]} // Evaluate, {t, t0 + .01, 100000 + .01}] /.
Solution // ReleaseHold // GraphicsRow


The Shooting Method is not strictly necessary here. In addition to MicheaelE2's solution:

t0 = 2;
M = 8;
f[t_] := t^2 - M/t;

Solution = NDSolve[
{x''[t]+(f'[t]/f[t]+2/t) x'[t]+(y[t]^2/f[t]^2+2/f[t])x[t] == 0,
y''[t] + 2/t y'[t] - 2 x[t]^2/f[t] y[t] == 0,
x[t0 + .01] == -3/2 t0 x'[t0 + .01],
y[t0 + .01] == 0,
y[100000 + .01] == 10,
x[100000 + .01] == 0},
{x, y}, {t, t0 + .01, 100000 + .01}, Method -> "ExplicitRungeKutta"];

Plot[Evaluate[{x[t], y[t]} /. Solution], {t, t0 + .01, 10},
PlotLegends -> {x[t], y[t]}]


The conditions:

-3/2 t0 x'[t0 + .01] /. Solution
x[100000 + .01] /. Solution


{1.03785} {-1.51313*10^-13}

y[t0 + .01] /. Solution
y[100000 + .01] /. Solution


{0.} {10.}

• This does the Shooting Method under the hood, which is the default for such BVPs. To specify the "StartingInitialConditions", which seemed to be the OP's intention, one has to do that through the Shooting Method specification. Commented Aug 17, 2015 at 16:57
• @Michael E2, that's right. The OP ask "... how should I solve it? I think the shooting method is required..." and his attempt to solve the DE was Method -> {"Shooting...."}. This Method definition in NDSolve is not strictly necessary.
– user31001
Commented Aug 17, 2015 at 17:25
• thanks @Michael E2 Willinski. Yeah it helps and solved the problem . But could you try to explain me the following thing... If I use shooting then I have to find out the best possible value of the derivative which makes the boundary value of the x[t] zero or the desired value. Whereas in the explicit Runga kutta method I don't need to put the value by hand. So does it have little more advantage in that sense or there is some other issues??
– kau
Commented Aug 17, 2015 at 17:30
• @kau I know the issues but I'm uncertain of the correct answer. Each time-integration method ("ExplicitRungeKutta", the default LSODA, etc.) has different step control, stiffness, error, etc. issues. Apparently the Runge-Kutta error etc. allows the shooting method to quickly find initial conditions that satisfy the BCs when integrating with that method, but not with LSODA. I doubt it always turns out that way. It must depend on the ODE. If LSODA is some other method is more accurate RK, you could use RK to get your "StartingInitialConditions"`. Commented Aug 17, 2015 at 18:39