# Finding the amplitude of function solved using ParametricNDSolve

I have been trying to find the peaks of a function I have plotted using ParametricNDSolve. I have to find these peaks to calculate the amplitude of all the various waves in the observed output. By amplitude here I mean the difference between an adjacent crest and trough. I shall attach the program here to clarify the problem.

Clear["Global*"]
om = 1;
k = 1;
L = 1/1000;
P = 1.3;

eqns = {a'[t] ==
I*((de/om)*a[t] - (b[t] + Conjugate[b[t]])*a[t] -
1/2) - (k/(2*om))*a[t],
b'[t] == -I*((P*(Abs[a[t]])^2)/2 + b[t]) - (L/(2*om))*b[t],
b[0] == 0, a[0] == 0};

s = ParametricNDSolve[eqns, {a, b}, {t, 0, 100}, de,
MaxSteps -> \[Infinity]];

x[de_][t_] = (b[de][t] + Conjugate[b[de][t]]) /. s;

Manipulate[
Plot[x[de][t], {t, 0, 99}, PlotRange -> {-4, 4}], {de, -3, 3, .1,
Appearance -> "Labeled"}]


Now I have to find the value of the amplitude of x for every value of de that I have plotted. While I have been able to find the peaks for individual values of de and a single peak in for these values, I have been unable to find all the peaks and amplitudes for all the values of de. I used various function including FindPeaks[] to try to find such values but was unable to do so. Normally I would have used for loops but that seems to advised against in Mathematica. I would also like to know whether such a code would be unwieldy and whether I use brute force to find such values.

I would be very grateful to anyone who could help me solve this dilemma. I would like to thank you for it beforehand

It might be easiest to detect the extrema within the NDSolve using WhenEvent.

Clear["Global*"]
om = 1; k = 1; L = 1/1000; P = 1.3;
tmax = 10000;

eqns = {
a'[t] == I*((de/om)*a[t] - (b[t] + Conjugate[b[t]])*a[t] - 1/2) - (k/(2*om))*a[t],
b'[t] == -I*((P*(Abs[a[t]])^2)/2 + b[t]) - (L/(2*om))*b[t],
b[0] == b0, a[0] == a0,
WhenEvent[b'[t] + Conjugate[b'[t]] == 0, Sow[{t, b[t] + Conjugate[b[t]]}]]
};

de = -3;
{a0, b0} = {0, 0};
{s, ex} = Reap[NDSolve[eqns, {a, b}, {t, 0, tmax}, MaxSteps -> \[Infinity]][[1]]];

(* plot saved extrema *)
ListPlot[ex[[1]], AxesLabel -> {"t", "extrema"}]


If you're interested in the long-term behavior, then it looks like tmax>100 is probably a good idea.

Then you can loop over de, saving the difference between the last two extrema as the amplitude (assuming a simple cycle) and using the final values of a and b as initial conditions for the next value of de.

Clear[de];
{a0, b0} = {a[tmax], b[tmax]} /. s;
tmax = 1000;
res = {};
Do[
{s, ex} = Reap[NDSolve[eqns, {a, b}, {t, 0, tmax}, MaxSteps -> \[Infinity]][[1]]];
{a0, b0} = {a[tmax], b[tmax]} /. s; (* save final values *)
AppendTo[res, {de, Abs[ex[[1, -1, 2]] - ex[[1, -2, 2]]]}];
, {de, -3, 3, 0.1}];

(* plot amplitude vs de *)
ListPlot[res, AxesLabel -> {"de", "amp"}]


That jump at de = 0.7 makes me think there might be alternative attractors. Indeed, if we generate the bifurcation diagram starting at de = 3 and moving towards lower de, we get a different result.

de = 3;
tmax = 10000;
{a0, b0} = {0, 0};
{s, ex} =  Reap[NDSolve[eqns, {a, b}, {t, 0, tmax}, MaxSteps -> \[Infinity]][[1]]];

Clear[de];
{a0, b0} = {a[tmax], b[tmax]} /. s;
tmax = 1000;
res = {};
Do[
{s, ex} = Reap[NDSolve[eqns, {a, b}, {t, 0, tmax}, MaxSteps -> \[Infinity]][[1]]];
{a0, b0} = {a[tmax], b[tmax]} /. s;
AppendTo[res, {de, Abs[ex[[1, -1, 2]] - ex[[1, -2, 2]]]}];
, {de, 3, -3, -0.1}];

ListPlot[res, AxesLabel -> {"de", "amp"}, PlotStyle -> Orange]


So a full understanding of the model's dynamics will require a bit of care!

pic=Plot3D[x[de][t], {t, 0, 99}, {de, -3, 3}, PlotRange -> {-4, 4},MaxRecursion -> 4]

If you split the region t,de NMaximize finds the global(!) peak in the subregion
NMaximize[{x[de][t], 0 < t < 10, -3 < de < 3}, {t, de}]