Mixing equations: numerical solutions

Notes: scale fields normalized to the total power first, then identify the characteristic length and use it to scale z

SHG

ξmax = 10;Δs=0.4;p1init=1;
ans = NDSolve[{
a1'[ξ]== I Conjugate[a1[ξ]] a2[ξ](E^(-I Δs ξ)),
a2'[ξ]==I ((a1[ξ])^2)(E^(I Δs ξ)),
a1[0]==(Sqrt[p1init]),
a2[0]==(Sqrt[1-p1init])},
{a1,a2},{ξ,0,ξmax}];
a1ans = a1 /.ans[[1,1]];
a2ans = a2 /.ans[[1,2]];

Plot[{((Abs[a1ans[ξ]])^2),((Abs[a2ans[ξ]])^2)},{ξ,0,ξmax},PlotStyle->{Red,Blue},PlotRange->{0,1}]