Starting with your equations we get
ode = y''[t] + y[t] + 0.25*y[t]^2 == 0.6*Sin[0.2*t];
sol = NDSolve[{ode, y'[0] == 0, y[0] == 1}, y, {t, 0, 300}];
Plot[Evaluate[y[t] /. sol], {t, 0, 250}, Frame -> True]
As mikado states this function does not appear to be decaying so we cannot take the full Fourier transform. However, we can take the Fourier transform of the section you have simulated. Working with this interval we can proceed to sample it at regular time steps.
sr = 3;(* Sample Rate*)
data = Table[Evaluate[y[t] /. sol[[1]]], {t, 0, 300, 1/sr}];
Here sr is the sample rate i.e. the number of samples per second. You need to choose this carefully. I have chosen 3 because you have about 50 cycles in 300 seconds so a frequency of about 0.17 Hz. It is necessary to have a sample rate larger than this to avoid aliasing and to look at higher frequencies than 0.17Hz (see Nyquist Frequency).
We now take the Fourier transform of the sampled data with an appropriate choice of FourierParameters .
ft = Fourier[data, FourierParameters -> {-1, -1}];
We also need a set of frequency values to provide abscissae for the spectral values which we calculate as follows.
ff = Table[(n - 1) sr/Length@data, {n, Length@data}];
Now we can plot the magnitude of the spectrum. I will leave you to calculate the phase.
ListLinePlot[Transpose[{ff, Abs[ft]}][[1 ;; 450]], PlotRange -> All,
Frame -> True, FrameLabel -> {"Frequency/Hz", "Spectral Level"}]
Fourier
. $\endgroup$