# Nonlinear differential equation system with initial conditions in the different point

I faced with unusual (at least for me) problems.

I have a three nonlinear equations

$$\dot x = f(x,y,z)\quad \dot y = g(x,y,z)\quad \dot z = h(x,y,z)$$

and the following initial conditions $$x(0)=1$$, $$y(0)=0$$, $$z(4)=1$$.

I try to solve it in the Mathematica

L = {D[Ap[z], z] == -gp/2 (Abs[Af[z]]^2 + Abs[Ab[z]]^2) Ap[z] +
I gpp (Abs[Ap[z]]^2 + 2 Abs[Af[z]]^2 + 2 Abs[Ab[z]]^2) Ap[z] -
alp/2 Ap[z],
D[Af[z], z] == -gs/2 Abs[Ap[z]]^2 Af[z] +
I gss (2 Abs[Ap[z]]^2 + Abs[Af[z]]^2 + 2 Abs[Ab[z]]^2) Af[z] +
I kk Ab[z] + I db Af[z] - a1s/2 Af[z], -D[Ab[z], z] ==
gs/2 Abs[Ap[z]]^2 Ab[z] +
I gss (2 Abs[Ap[z]]^2 + 2 Abs[Af[z]]^2 + Abs[Ab[z]]^2) Ab[z] +
I kk Af[z] + I db Ab[z] - a1s/2 Ab[z]}

vp = 1.; vs = 2; gp = 3; gs = 4; kk = 5; a1s = 0.5; gpp = 0.1; gss = 0.23; alp = 0.79; db = 0.07;

s = NDSolve[Join[L, {Ap[0] == 1, Af[0] == 0, Ab[4] == 1}], {Ap, Af, Ab}, {z, 0, 4}]


But the solution $$s$$ is not satisfied the condition $$Ab[4] =1$$. In the Mathematica solution $$Ab[4]=-7.27787 + 1.55824 I$$ instead of $$1$$.

How to solve it correctly?

• Is the solution defined in the complex plane Z or are only complex functions Ap, Af, Ab of real variable z? Commented Oct 7, 2018 at 5:38
• $z$ is a real. Functions $Ap, Af, Ab$ are complex. Commented Oct 7, 2018 at 5:47

First, we reduce the system of equations to a real form. Secondly, we will define the initial data using a parametric solution.

L = {-D[Ap[z], z] - gp/2 (Abs[Af[z]]^2 + Abs[Ab[z]]^2) Ap[z] +
I gpp (Abs[Ap[z]]^2 + 2 Abs[Af[z]]^2 + 2 Abs[Ab[z]]^2) Ap[z] -
alp/2 Ap[z], -D[Af[z], z] - gs/2 Abs[Ap[z]]^2 Af[z] +
I gss (2 Abs[Ap[z]]^2 + Abs[Af[z]]^2 + 2 Abs[Ab[z]]^2) Af[z] +
I kk Ab[z] + I db Af[z] - a1s/2 Af[z],
D[Ab[z], z] + gs/2 Abs[Ap[z]]^2 Ab[z] +
I gss (2 Abs[Ap[z]]^2 + 2 Abs[Af[z]]^2 + Abs[Ab[z]]^2) Ab[z] +
I kk Af[z] + I db Ab[z] - a1s/2 Ab[z]};

Ap'[z] = Ap1'[z] + I*Ap2'[z];

Af'[z] = Af1'[z] + I*Af2'[z];
Ab'[z] = Ab1'[z] + I*Ab2'[z];
Ap[z] = Ap1[z] + I*Ap2[z];
Af[z] = Af1[z] + I*Af2[z];
Ab[z] = Ab1[z] + I*Ab2[z];

L // FullSimplify;
eq = ComplexExpand[L];
eqs =
Table[{ComplexExpand[Re[eq[[i]]]] == 0,
ComplexExpand[Im[eq[[i]]]] == 0}, {i, 1, 3}];

vp = 1.; vs = 2; gp = 3; gs = 4; kk = 5; a1s = 0.5; gpp = \
0.1; gss = 0.23; alp = 0.79; db = 0.07;

AB =
ParametricNDSolveValue[
Join[eqs, {Ap1[0] == 1, Af1[0] == 0, Ab1[0] == p, Ap2[0] == 0,
Af2[0] == 0, Ab2[0] ==q}], Ab1, {z, 0, 4}, {p,q}];
AC = ParametricNDSolveValue[Join[eqs, {Ap1[0] == 1, Af1[0] == 0, Ab1[0] ==
p, Ap2[0] == 0, Af2[0] == 0, Ab2[0] == q}], Ab2, {z, 0, 4}, {p, q}];
FindRoot[{AB[p, q][4] == 1, AC[p, q][4] == 0}, {p, -0.5,
0}, {q, -.5, 0}, Method -> "Secant"]

Out[]= {p -> 4.3433*10^-8, q -> 5.71603*10^-9}

{Plot[Evaluate[{AB[4.3433002424791235*^-8, 5.716028592104211*^-9][
z], AC[4.3433002424791235*^-8, 5.716028592104211*^-9][z]}], {z,
0, 4}, PlotRange -> All],
Plot[Evaluate[{AB[4.3433002424791235*^-8, 5.716028592104211*^-9][
z], AC[4.3433002424791235*^-8, 5.716028592104211*^-9][z]}], {z,
0, 4}, PlotRange -> Automatic, PlotLegends -> Automatic]}


Finally, we find the entire solution.

s = NDSolve[
Join[eqs, {Ap1[0] == 1, Af1[0] == 0,
Ab1[0] == 4.3433002424791235*^-8, Ap2[0] == 0, Af2[0] == 0,
Ab2[0] == 5.716028592104211*^-9}], {Ap1, Ap2, Af1, Af2, Ab1,
Ab2}, {z, 0, 4}]
{Plot[Evaluate[{Ap1[z], Ap2[z], Af1[z], Af2[z], Ab1[z], Ab2[z]} /.
s], {z, 0, 4}, PlotRange -> {-.1, 0.3}],
Plot[Evaluate[{Ap1[z], Ap2[z], Af1[z], Af2[z], Ab1[z], Ab2[z]} /.
s], {z, 0, 4}, PlotLegends -> Automatic, PlotRange -> All]}


• Thank you! I understand the idea. But realization of this idea is not right. Because in the function ParametricNDSolveValue should be to parameters $p1$ and $p2$. Ab1[0]=p1 and Ab2[0]=p2. In the functions FindRoot will be two equation $Ab1[4]==1$ and $Ab2[4]==0$. But I faced with the problem that FindRoot does not give me a result and says "FindRoot::nlnum: "The function value {-1.+0.1[4.],0.2[4.]} is not a list of numbers with dimensions {2} at {p1,p2} = {0.1,0.2}."" Commented Oct 7, 2018 at 15:20
• @Peter In your problem, there are no boundary conditions for the imaginary parts of the functions. Think about what the conditions should be, and then I will help solve the problem. Commented Oct 7, 2018 at 15:42
• In my problem Ab(4)=1. It means that Ab1(4)=1 and Ab2(4)=0. In your answers Ab2(4) has non zero value. Commented Oct 7, 2018 at 15:52
• Are other boundary conditions correct? Commented Oct 7, 2018 at 16:12
• Yes. The boundary conditions are $Ap[0] == 1, \quad Af[0] == 0,\quad Ab[4] == 1$. In the your notations are $Ap1[0] == 1,\quad Af1[0] == 0, \quad Ab1[4] == 1 \quad Ap2[0] == 0, \quad Af2[0] == 0, \quad Ab2[4] == 0$ Commented Oct 7, 2018 at 16:37