fizyka XD

Clear[x]
indeks = 145406
f = 1
\[Omega]0 = indeks + 1
b = 1/4 * \[Omega]0
\[Omega]1 =(\[Omega]0 ^2 - 1/2 * b^2)^(1/2)
\[Omega]2 = 3/4 * (\[Omega]0^2 - 1/2 * b^2)^(1/2)
\[Omega]3 = 5/4 * (\[Omega]0^2 - 1/2 * b^2)^(1/2)
s[\[Omega]_] := NDSolve[{x''[t] + b x'[t] + \[Omega]0^2 x[t] == f Sin[\[Omega] t], x[0] == 0, x'[0] == 0}, x[t], {t, 0, 2*Pi/ \[Omega]0}]
x[t_, \[Omega]_] := s[\[Omega]][[1, 1, 2]]

Plot[Evaluate[{x[t, \[Omega]1]}], {t, 0, 2*Pi/ \[Omega]0}, PlotRange -> Automatic,PlotStyle ->Automatic, AxesLabel -> {"t", "x[t]"}]

Plot[Evaluate[{x[t, \[Omega]2]}], {t, 0, 2*Pi/ \[Omega]0}, PlotRange -> Automatic,PlotStyle ->Automatic, AxesLabel -> {"t", "x[t]"}]

Plot[Evaluate[{x[t, \[Omega]3]}], {t, 0, 2*Pi/ \[Omega]0}, PlotRange -> Automatic,PlotStyle ->Automatic, AxesLabel -> {"t", "x[t]"}]

xs = Range[19];
ys = Range[19];
For[k = 1, k <= 19, k++,
xs[[k]] = k/10 * (\[Omega]0^2 - 1/2 * b^2)^(1/2);
ys[[k]] = First[NMaximize[{x[t, k/10 * (\[Omega]0^2 - 1/2 * b^2)^(1/2)], 0 <= t <= 2*Pi/\[Omega]0}, {t}]]
]
ListPlot[Thread[{xs,ys}]]

m= Table[i,{i  ,{xs,ys}}];
Grid[Transpose@m,Frame ->All]

xmax = Max[ys]    (* wartosc X0 *);
x0\[Omega]max = Max[xs]  (* wartosc omega max *);
\[Omega]minus = xs[[1]](*wartosc omega minus*);
\[Omega]plus = xs[[20]] (*wartosc omega plus*);
\[CapitalDelta]\[Omega] = \[Omega]plus - \[Omega]minus (*wartość omega delta*);

N[\[Omega]minus]
N[\[Omega]plus]
N[\[CapitalDelta]\[Omega]]