Solution to differential equation only starts at random times (numerical error?)

I am solving a set of differential equations with an oscillatory forcing (square wave). Instead of the solution starting to oscillate immediately it stays still until some time. That varies with small changes to initial conditions or parameters. An image will exemplify better. It also has this weird behavior if I change the initial angle -pi/2 by 0.1 for example.

Has anyone encountered this? It seems impossible to me that this comes from the equations, Debugging suggestions would be useful.

The code:

totaltime = 20;
\[Phi]e[t_] := SquareWave[{phi0, Pi/2}, t*freq];
efieldm[t_] := SquareWave[{m0, 1}, t*freq] ;
eqphi := \[Gamma] Derivative[\[Phi]bg][t] ==
1/2 (-efieldm[t]^2) Sin[2 (\[Phi]bg[t] - \[Phi]e[t])];

efOg = {1, 3};
phi0 = ArcTan[efOg[]/efOg[]];
m0 := Norm[efOg]
freq = 0.5;
v2 := NDSolve[
eqphi && \[Phi]bg == -(\[Pi]/2), \[Phi]bg[t], {t, 0,
totaltime}];
start1 = 0.2; n1 = 5; step1 = 0.15;
table1 = Table[v2, {freq, start1, start1 + n1*step1, step1}];
table2 = Table[efieldm[t], {freq, start1, start1 + n1*step1, step1}]
For[n = 1, n <= 5, n++,
Print[Plot[{\[Phi]bg[t] /. table1[[n]], table2[[n]]}, {t, 0,
totaltime}]]]

• \[Gamma] seems undefined. I can't run the code. Nov 18 '21 at 15:39
• Try limiting step size. E.g. MaxStepSize -> 0.1. Nov 18 '21 at 15:41
• Or try PiecewiseExpand@SquareWave[..], which will create events for the discontinuities. Basically, the ode behaves so well at first that NDSolve takes too big a step that bypasses some oscillations of SquareWave. (Consider reporting it to WRI. NDSolve should behave better.) Nov 18 '21 at 15:42
• I vote for the bugs tag: "NDSolve automatically does processing for discontinuous functions like Sign" (from the docs), and SquareWave is like Sign, imo. Nov 18 '21 at 15:49
• @MichaelE2, Sorry about gamma, its just a parameter. And excelent sugestion! It corrected all the variations I tried. I was already trying python, but now this works. Nov 18 '21 at 15:53