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]]
x[t, \[Omega]1]

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