Trying to solve non-linear coupled differential equation with boundary conditions at different points

Question

I'm trying to solve systems of dif. equations with boundaries conditions at distinct points such as:

inf = 100;
test1 = x''[t] - x[t] + 3 x[t] y[t] == 0;
test2 = y''[t] - y[t] - 2 (x[t] + y[t])^2 == 0;

sol = {x[t], y[t]} /. 
  First[NDSolve[{test1, test2, x[inf] == 2, y[inf] == 4, x[-inf] == 1,
      y[-inf] == 2}, {x[t], y[t]}, {t, -inf, inf}]]

I get no solution for this set of boundary conditions. The curious thing is that with a simple change (ex. y[inf] == 0) the program returns me the IntepolatingFuntion solutions as expected. I'm trying to find the reason for this behavior. I wonder if NDSolve is the apropriate function to deal with a problem with these boundary constraints. I hope that someone can point me to a more fruitful approach to the problem or help me understand why this is happening.

Answer 1

About the largest value of inf that yields a result for the code as given in the question is 0.6.

SetOptions[Plot, ImageSize -> Large, LabelStyle -> {15, Bold, Black}];

inf = 3/5; wp = 30;
sol = NDSolveValue[{test1, test2, x[inf] == 2, y[inf] == 4, 
    x[-inf] == 1, y[-inf] == 2}, {x[t], y[t]}, {t, -inf, inf}, 
    WorkingPrecision -> wp];
g = D[sol, t] /. t -> -inf
Plot[sol, {t, -inf, inf}]

(* {2.33385632500390682864115515437, -6.75800477617055150858467720140} *)

To increase inf, it is necessary to use the "Shooting" Method explicitly. Estimating good initial guesses a priori is, however, difficult. One approach is to increase inf gradually, using the results from a solution for inf to estimate the shooting initial guesses for a slightly larger value of inf. Progressively larger WorkingPrecision also is needed. For instance,

Do[sol = 
   NDSolveValue[{test1, test2, x[inf] == 2, y[inf] == 4, x[-inf] == 1, y[-inf] == 2}, 
      {x[t], y[t]}, {t, -inf, inf}, WorkingPrecision -> wp, 
      Method -> {"Shooting", "StartingInitialConditions" -> 
      Thread[{x'[-inf], y'[-inf]} == g]}]; 
   g = SetPrecision[D[sol, t] /. t -> -inf, wp]
    , {inf, 6/10, 8, 1/10}]
g
Plot[sol, {t, -pinf, pinf}, PlotRange -> All, AxesLabel -> {t, "x,y"}]
sol /. t -> pinf

(* {0.425811455218237081374463614922, -6.11260768006832031438841518769} *)

(* {1.8942478267433165614815184445, 3.976808432427045469940175848} *)

As can be seen from sol /. t -> pinf, the accuracy with which the pinf boundary conditions are satisfied is deteriorating for inf = 15. There are, of course, ways to improve the accuracy, particularly by using a larger WorkingPrecision. Nonetheless, as inf increase further and the dependent variables vary ever more rapidly near the boundaries, the computation will become ever slower. The simple approach shown here provides a picture of what is happening to the solution and how to extend it to moderately large inf.