Just use ParametricNDSolveValue to solve your problem:
sol = ParametricNDSolveValue[{x''[t] +2 \[Delta]*x'[t] + \[Omega]0^2*x[t] +\[Alpha] x[t]^2 == 0,x[0] == 1, x'[0] == 0}, x, {t, t0, t1}, {\[Alpha]}]
(*ParametricFunction[ <> ] *)
ParametricPlot[Evaluate[Table[{sol[\[Alpha]][t],
sol[\[Alpha]]'[t]}, {\[Alpha], {\[Alpha]1, \[Alpha]2, \[Alpha]3}}]], {t, t0,t1}, PlotRange -> All]
