# Inaccuracy with this solution of this PDE from Mathematica NDSolve

I have a problem with this script that solves a nonlinear PDE. Please note that the PDE is defined for the complex function \Psi of x and t whose module is called \phi and argument called \chi.

ClearAll["Global*"];
M = 1; (* Mass scale*)
A = M*10^-3; (* Gaussian Amplitude*)
L = 10/M; (* Length scale*)
Lb = 100*L; (* Boundary distance*)
tmax = 30*L; (* Maximal time*)
\[Chi]0[x_] := A*Exp[-x^2/(2*L^2)]; (* IC for light field*)
\[Phi]0[x_] := Sqrt[M^2 - 0 + (\[Chi]0'[x])^2/M^2]; (* IC for heavy field*)
sol = NDSolve[{
-D[\[CapitalPsi][x, t], t, t] + D[\[CapitalPsi][x, t], x, x] +
2*((Abs[\[CapitalPsi][x, t]])^2 - M^2 ) \[CapitalPsi][x, t] == 0,
\[CapitalPsi][x, 0] == \[Phi]0[x]*Exp[I*\[Chi]0[x]],
(D[\[CapitalPsi][x, t], t] /. t -> 0) == 0,
D[\[CapitalPsi][-Lb, t], t] == D[\[CapitalPsi][Lb, t], t],
\[CapitalPsi][-Lb, t] ==  \[CapitalPsi][Lb, t]
},
{\[CapitalPsi]}, {x, -Lb, Lb}, {t, 0, tmax}];

y[x_, t_] = \[CapitalPsi][x, t] /. sol[];
\[Chi][x_, t_] = Arg[y[x, t]];
\[Phi][x_, t_] = Abs[y[x, t]];
Plot[{\[Chi]0[x], \[Chi][x, 0]}, {x, -Lb, Lb}, PlotRange -> All]


The last line helps me to compare the function that defines the initial condition for \chi to the evaluation of the evaluation of the argument of \Psi at time t=0 (these two are supposed to match). Unfortunately, the plot gives a different result for the two functions (the real initial condition in blue and the solution in orange): I don't understand why the solution does not respect the initial condition. Can anybody help me please?

Many thanks.

It looks like the solver has trouble converging. If you use smaller step sizes or other measures to increase numerical stability it should help. For example with MaxStepSize -> 0.5:
ClearAll["Global*"]; M = 1;(*Mass scale*)A = M*10^-3;(*Gaussian Amplitude*)L = 10/M;(*Length scale*)Lb = 100*L;(*Boundary distance*)tmax = 30*L;(*Maximal time*)\[Chi]0[x_] := A*Exp[-x^2/(2*L^2)];(*IC for light field*)\[Phi]0[x_] := Sqrt[M^2 - 0 + (\[Chi]0'[x])^2/M^2];(*IC for heavy field*)sol = NDSolve[{-D[\[CapitalPsi][x, t], t, t] + D[\[CapitalPsi][x, t], x, x] + 2*((Abs[\[CapitalPsi][x, t]])^2 - M^2) \[CapitalPsi][x, t] == 0, \[CapitalPsi][x, 0] == \[Phi]0[x]* Exp[I*\[Chi]0[x]], (D[\[CapitalPsi][x, t], t] /. t -> 0) == 0, D[\[CapitalPsi][-Lb, t], t] == D[\[CapitalPsi][Lb, t], t], \[CapitalPsi][-Lb, t] == \[CapitalPsi][Lb, t]}, {\[CapitalPsi]}, {x, -Lb, Lb}, {t, 0, tmax}, MaxStepSize -> 0.5];
y[x_, t_] = \[CapitalPsi][x, t] /. sol[]; \[Chi][x_, t_] = Arg[y[x, t]]; \[Phi][x_, t_] = Abs[y[x, t]]; Plot[{\[Chi]0[x], \[Chi][x, 0]}, {x, -Lb, Lb}, PlotRange -> All]` 