Nonlinear differential equation system with initial conditions in the different point

Question

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?

Answer 1

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]}