Just use NDSolveValue 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]
[![enter image description here][1]][1]
    


  [1]: https://i.sstatic.net/RxfQh.jpg