# How to Plot the Steady State Phase Response Curve?

I have a code to plot the steady state frequency response curve (provided by @ Michael E2), How can I plot the steady state phase response with the help of this code same as shown in the figure:

.

\[Kappa]=1.4;
WE = 0.2;
\[Omega] = Sqrt[3*\[Kappa]*(1 + WE) - WE];
f = 0.2;
\[Epsilon] = 0.2;

amp[w_, w0_?NumericQ, damping_?NumericQ] /; w == 0 = 1.;
amp[w_?NumericQ, w0_?NumericQ, damping_?NumericQ] :=
Block[{x}, #[#["Domain"][[1, -1]]] &[
x /. First@
NDSolve[{x''[
T] == (1/\[Epsilon])*((((1 +
WE)*(1 + \[Epsilon]*x[T])^(-3*\[Kappa] -
1))) - ((1 + \[Epsilon]*
x[T])^-1) - (WE*(1 + \[Epsilon]*x[T])^-2) - (4*
damping*\[Epsilon]^2*(x'[
T])*(1 + (\[Epsilon]*
x[T]))^-2) - ((3*\[Epsilon]^2*(x'[
T]^2)*(1 + \[Epsilon]*x[T])^-1)/2) - \[Epsilon]^3*
f*Sin[T* (w0 + \[Epsilon]^2 w)]), x[0] == 1, x'[0] == 0,
WhenEvent[x'[T] < 0 && T > 5000, "StopIntegration"]},
x, {T, 0, 5010}]]];
Off[NDSolve::ndsz]
pltNum = Plot[
Evaluate[Table[
amp[w, \[Omega], z], {z, {\[Epsilon]*0.04}}]], {w, -5, 5},
PlotStyle -> {Thick, Dashed, Red},
PlotLegends -> Placed[{"Numerical"}, {0.8, 0.8}], PlotRange -> All]


Any help would be highly appreciated. Thanks!