# NDSolve with one initial condition and one terminal condition

I have been trying to solve two coupled first-order differential equations. For one of those two equations, the initial condition is known and for the other one, the terminal condition is known. However, Mathematica is showing error for such type of mixed conditions. For example, consider the following simple example -

xi = 0;
xf = 1;
sol1 = NDSolve[{G'[x] == 0, H'[x] == 0, G[xi] == 0, H[xf] == 1}, {G[x], H[x]}, {x, xi, xf}]


The error message is -

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0..


Is there a way to solve this issue?

• Maybe implement a shooting method? Commented Oct 16, 2022 at 15:34
• You still need an initial condition for H[x] ! Commented Oct 16, 2022 at 15:37
• G'[x]==0 means G is a constant. Similar for H. Therefore, G[x]==0 and H[x]==1. Commented Oct 16, 2022 at 15:46
• @PoreyS I've added an answer with an example. The ODEs you posted where too trivial, so I changed the system, but it could be generalized to what you want. Commented Oct 16, 2022 at 18:32
• @UlrichNeumann The OP has initial conditions, albeit at different points and the goal is to solve the ODEs with those. Commented Oct 16, 2022 at 18:35

As the ODEs posted by the OP are maybe too trivial for this, I took the liberty and changed them to a really coupled system (taken straight from the Help) in order to highlight the procedure.

The idea of the "shooting method" I mentioned in the comments is to try to find the initial condition for the 2nd variable, in this example $$y(t)$$, such that the requested terminal condition is satisfied.

So, let's assume the following system with initial conditions at $$t_i=0$$ $$x(t_i)=1$$ and $$y(t_i)=y_0$$, while the terminal condition at $$t_f=20$$ is $$y(t_f)=0.15$$. Then we set up the solution with NDSolve as a free parameter for $$y_0$$ as:

ti = 0;
tf = 20;
yf = 0.15;
sol[y0_?NumericQ] := NDSolve[{x'[t] == -y[t] - x[t]^2, y'[t] == 2 x[t] - y[t]^3, x[ti] == 1, y[ti] == y0}, {x, y}, {t, ti, tf}]


Let us now define the value of $$y(t)$$ at $$t_f$$ as a function of the initial condition $$y_0$$ and use FindRoot to solve for $$y_0$$:

Yf[y0_?NumericQ] := y[t] /. sol[y0][[1]] /. t -> tf
ic = FindRoot[Yf[y0] == yf, {y0, -1, -0.9}, MaxIterations -> 1000]
(* {y0 -> -0.612372} *)


Indeed this confirms that the value gives the proper "terminal condition":

Yf[ic[[1, 2]]]
(* 0.15 *)


Also, plotting the solution proves the same (the red point is the "terminal condition"):

soly0 = sol[ic[[1, 2]]][[1]];
ParametricPlot[{x[t], y[t]} /. soly0, {t, ti, tf}, Frame -> True, FrameLabel -> {"x(t)", "y(t)"}, FrameStyle -> Black,Epilog -> {Red, Point[{x[t], y[t]} /. soly0 /. t -> tf]},PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}, PlotStyle -> Black]


However, caution is needed as many times there may be more than one initial conditions that give the same terminal condition as shown by the plot of $$y_f$$ vs $$y_0$$, as the function $$y_f(y_0)$$ in general may not be monotonic (as in this example):

Intuitively, this is due to the fact that more than one trajectories may pass from the same point, i.e. the terminal condition, for different initial conditions.

• Thanks for the reply. I have two questions - 1. Yf[y0_?NumericQ] := y[t] /. sol[y0][[1]]... what does sol[y0][[1]] stands for? 2. {y0, -1, -0.9} -> what does this range for? Commented Oct 17, 2022 at 2:54
• The first one is the numerical solution of NDSolve for each value of the initial condition, while the latter is just a starting range to help FindRoot. Commented Oct 17, 2022 at 6:21
• @PoreyS If you think my answer addresses your problem, please consider accepting it or clarifying what else you need, eg maybe an ODE system closer to what you want. Commented Oct 17, 2022 at 8:36
• Thanks a lot. I was trying to use your code for my scenario. Commented Oct 17, 2022 at 18:18

First of all, these are not coupled. If you solve them separately you will get correct solutions.

NDSolve got confused since it wanted to know what H value at zero and you did not tell. So it complained. It did not know these could be solved separately. After all, it is numerical solver.

Clear["Global*"]
xi = 0;
xf = 1;
solG = NDSolve[{G'[x] == 0, G[xi] == 0}, G, {x, xi, xf}]
solH = NDSolve[{H'[x] == 0, H[xf] == 1}, H, {x, xi, xf}]
p1 = Plot[Evaluate[G[x] /. solG], {x, xi, xf}];
p2 = Plot[Evaluate[H[x] /. solH], {x, xi, xf}];
Show[p1, p2]


These solutions match the analytical DSolve solution. I do not know enough about NDSolve to know if there is an option to avoid this. May be there is.

DSolve[{G'[x] == 0, H'[x] == 0, G[xi] == 0, H[xf] == 1}, {G[x],H[x]}, x]


So the rule of the thumb, is that if you are solving numerically more than one ode at same time (coupled or not), they all have to have same initial conditions (default behavior).

Otherwise, different setting/options might need to be added explicitly.

Comment says to try shooting method. But I do not know how to do this with two ode's one with initial conditions and the other without that. I do not know even if this is supported. Someone else would know for sure.