# Jumps in NDSolve results

I need to compute using NDSolve routine, some function $F(x)$, having two possible values $F_1(x)$ and $F_2(x)$ depending on whether the argument exceeds some critical value $x>x_c$. The problem is that NDSolve routine return "jumping results" for the region $x>x_c$.

As an example of a problem assume following system (some simple forced pendulum with nonlinear frequncy shift, so-called "fold over" effect):

 w00 = 7000;
G0 = 50;
Q = 14000;
hU = 0.6*G0;
w0=w00-x

e1:=D[A[t], t] + G0*A[t] + I*(w0+Q*Abs[A[t]]^2 )*A[t] ==I*hU;

NDSolve[{e1, A[0] ==10^-6, {A[t]}, {t, 0, tmax},
MaxSteps -> Infinity, AccuracyGoal -> 50, MaxStepSize -> 0.01];


I am interested in $|A(t)|^2$ dynamics depending on the parameter x. Here $G0,Q,hU$ are some constants

Exact solution:

    p4 = ContourPlot[{A == (hU^2/((G0)^2 + (w0 - x + Q*A)^2))}, {x, 6000,
9000}, {A, 0, 0.2}, PlotPoints -> 40, ContourStyle -> {Dashed, Thick}];
Show[p4,PlotRange -> All]

Numeric solution using module:

Calc[x_] := Module[{q = x},
w0 = 7000;
G0 = 50;
Q = 14000;
hU = 0.6*G0;

w00 = w0 - q + Q*Abs[A0[t]]^2;
At = 10^-6;
G1 = 1.0 G0;
tmax = 4;
e1 := D[A0[t], t] + G0*A0[t] + I*w00*A0[t]== I*hU;

s2 = NDSolve[{e1,
A0[0] == (0.7 + 0.5 I)  + 0*10^-6}, {A0[t]}, {t, 0, tmax},
MaxSteps -> Infinity, AccuracyGoal -> 50, MaxStepSize -> 0.01];

{Evaluate[Abs[A0[t] /. s2 /. t -> tmax]][[1]]
}

]
h1 = {};
Dynamic[z]
For[z = 9000, z > 6000, z = z - 50,
temp = Calc[z];
h1 = Append[h1, {z, temp[[1]]^2}];
]
num = ListPlot[{h1}, PlotRange -> All, Joined -> True,
PlotStyle -> Thick]
Show[num]


Comparing both:

Show[p4,num]


Solid line – NDSolve solution, dashed line - exact solution

The result and number of jumps greatly depends on initial conditions.

A[0]=0.7+0.5*I


A[0]=10^-6


The problem is I want to divide the two "branches" of the function $F_1$ and $F_2$ and plot them independently. How can I do this using NDSolve and force it to go along the chosen branch and not "jump?"?

-