In my opinion the undamped case for a driven oscillator is a bit more subtle, because for the damped one you can argue that after a sufficiently long time (long against the inverse of the friction coefficient) after switching on the harmonic force the oscillator is in a stationary state oscillating with the frequency of the imposed harmonic force.
In the undamped case in general you have a superposition of the undamped natural oscillations and the ones imposed by the external harmonic force. Of course, as already mentioned above, there is the case where the frequency of the harmonic force coincides with the eigenfrequency of the oscillator, and then the amplitude grows linearly with time ("resonance catastrophe").
That said, let's do the calculation for the undamped oscillator. It's a bit simpler using the complex representation, i.e., introduce a complex variable ##z(t)## with the meaning that the physical quantity under consideration is ##x(t)=\text{Re} z(t)##. Then we want to solve the equation of motion
$$\ddot{x}+\omega_0^2 x=A \exp(\mathrm{i} \omega t).$$
The general solution is the general solution of the homogeneous equation (with the rhs of the equation set to 0):
$$z_{\text{hom}}(t)=C_1 \cos(\omega_0 t) + C_2 \sin(\omega_0 t),$$
and an arbitrary special solution of the inhomogeneous equation.
It's intuitive to make the ansatz
$$z(t)=B \exp(\mathrm{i} \omega t).$$
As we'll see in a moment that always works, except for the resonance case. So let's first assume ##\omega \neq \omega_0##. Plugging the ansatz into the equation we get
$$-(\omega^2-\omega_0^2) B = A \; \Rightarrow \; B=\frac{A}{\omega_0^2-\omega^2}.$$
Since both ##A## and ##B## are real, there's no phase shift between the external force and and that special solution of the inhomogeneous equation. I.e., you get
$$x_{\text{inh}}=\mathrm{Re} [B \exp(-\mathrm{i} \omega t)]=B \cos(\omega t) = \frac{A}{\omega_0^2-\omega^2} \cos(\omega t).$$
For the general solution
$$x(t)=C_1 \cos(\omega_0 t) + C_2 \sin(\omega_0 t) +\frac{A}{\omega_0^2-\omega^2} \cos(\omega t),$$
it doesn't make sense to talk about a phase shift, because it's a superposition of two harmonic motions with different frequencies. Only if ##C_1=C_2=0## you get a harmonic motion with the driving frequency ##\omega## and then there's no phase shift between ##x## and the external force.
The resonance case needs extra treatment. Here the idea is to use the ansatz of the rhs. of the equation with variation of the constant, i.e.,
$$z(t)=B(t) \exp(\mathrm{i} \omega_0 t).$$
Using
$$\ddot{z}=(\ddot{B} + 2 \mathrm{i} \omega_0 \dot{B} -\omega_0^2 B) \exp(\mathrm{i} \omega t),$$
in the equation of motion we find
$$\ddot{B} + 2 \mathrm{i} \omega_0 \dot{B} =A.$$
Since we only need one particular solution, we can solve this by making the ansatz
$$B(t)= B_0 t$$
leading to
$$B_0=\frac{A}{2 \mathrm{i} \omega_0}.$$
Thus the physical particular solution is
$$x_{\text{inh}}(t)=\mathrm{Re} B_0 t \exp(\mathrm{i} \omega_0 t) = \frac{A}{2 \omega_0} t \sin(\omega_0 t).$$
So the amplitude grows linearly with ##t##, and for large times you can neglect the solutions of the homogeneous equations. The phase shift of the oscillations is obviously ##-\pi/2##, because ##\sin(\omega_0 t)=\cos(\omega_0 t-\pi/2)##, i.e., the phase of the oscillation of ##x## is by 90 degrees behind the phase of the external force, but still there's the growing amplitude of this oscillation!
