Forced oscillation resonance frequencies Given this forced oscillator:
$$
\ddot{x}+\gamma\dot{x}+\omega_0^2x=\frac{F(t)}{m}
$$
Where $F(t)$ is:
$$
F(t)=\sum_{n=1}^{\infty}\frac{4F_0}{n\pi}\sin\left(\frac{2n\pi t}{T}\right)
\hspace{0.4cm} \text{  for }
n = 1, 3, 5, ...
$$
I start with
$
x_s(t)= A\sin(\omega t)+B\cos(\omega t)
$
And I solve them for A and B (see below), if I am not mistaken:
$$
A=\frac{F(\frac{\pi}{2\omega})}{m}
\frac{\left( \omega^2-\omega^2_0\right)}{\gamma^2\omega^2+\left( \omega^2-\omega^2_0\right)}
\hspace{0.2cm},\hspace{0.2cm}
B=\frac{F(\frac{\pi}{2\omega})}{m}
\frac{\gamma\omega}{\gamma^2\omega^2+\left( \omega^2-\omega^2_0\right)}
$$
Is this solution valid?

How can I calculate the resonance frequency? Because it seems that it is not $\omega=\omega_0$

Is there more than one resonance frequency?




My calculations
$$
\dot{x}_s=\omega A\cos(\omega t)-\omega B \sin(\omega t)
~~,~~~~
\ddot{x}_s=-\omega^2 A\sin(\omega t)-\omega^2 B \cos(\omega t)
$$
$$
\ddot{x}+\gamma\dot{x}+\omega_0^2x=\frac{F(t)}{m}\\
\!\!\!\! 
=-\omega^2 A\sin(\omega t)-\omega^2 B \cos(\omega t)
+\gamma\omega A\cos(\omega t)-\gamma \omega B \sin(\omega t)
+\omega_0^2A\sin(\omega t)+\omega_0^2B\cos(\omega t),
$$
and I solve for A and B  at $t=0$ and $t=\frac{\pi}{2\omega}$.
 A: Absorb the superfluous parameter m into your r.h.s. driving term, an obvious square wave with an infinity of harmonics,
$$
\ddot{x}+\gamma\dot{x}+\omega_0^2x=  \sum_{n=1,3,5...}^{\infty}\frac{4F_0}{n\pi}\sin\left( n\omega t\right),$$ where I've defined $\omega\equiv 2\pi/T$,
not your Ansatz parameter.
From the linearity of the equation, and the ultimate vanishing of the  transient solution of the damped undriven component, you realize you only need solve one problem, e.g. the leading harmonic ω, replicate/adapt the solution for each one nω, and add all of them together (an infinity) to get the most general solution of the full problem.
So you really only need solve the customary
$$\ddot{x}+\gamma\dot{x}+\omega_0^2x=   \frac{4F_0}{\pi}\sin\left( \omega t\right),$$
which you started doing with an ansatz of frequency ω, except for the full driving term, which cannot work. The steady-state solution is, then, taking the t=0 limit into consideration,
$$
x= -\frac{4F_0}{\pi}\left ( \sin(\omega t)\frac{\omega^2-\omega_0^2}{(\omega^2-\omega_0^2)^2+\gamma^2\omega^2}+ \cos(\omega t)\frac{\gamma\omega}{(\omega^2-\omega_0^2)^2+\gamma^2\omega^2} \right )\\ \equiv
\frac{4F_0}{\pi\sqrt{(\omega^2-\omega_0^2)^2+\gamma^2\omega^2}}\bigl ( \sin(\omega t)\cos \phi  + \cos(\omega t)\sin \phi   \bigr ) \\
= \frac{4F_0}{\pi\sqrt{(\omega^2-\omega_0^2)^2+\gamma^2\omega^2}} \sin(\omega t + \phi)   ~  , \leadsto \\
\bullet ~~~~~~~~\tan \phi ={ \gamma \omega\over \sqrt{(\omega^2-\omega_0^2)^2+\gamma^2\omega^2}. }~, 
$$
as customary.
The minimum of the square root denominator of the amplitude is at
$$
\omega=\sqrt{\omega_0^2-\gamma^2/2}~,
$$
the resonant frequency of a single harmonic forcer. It is slightly below the natural frequency $\omega_0$ for weak friction.
Now repeat all this (in your mind) for a single forcing term of the form
$ \frac{4F_0}{n\pi}\sin\left(n \omega t\right)$, and add all these solutions to the above, by linearity, to get the general steady-state solution (that is, ignore the transient). Note each n will have a characteristic, different phase, and a characteristic resonance frequency at 1/n of the above resonant frequency. Each term will be suppressed by 1/n. Because of the different phases, all near π/4 at resonance, there won't be any mode-locking making these lower resonant frequencies invisible. In practice, you see the fundamental harmonic resonance slightly below  $\omega_0$, but also the next odd harmonic resonance at somewhat below  $\omega_0/3$ !
I'm sure there are better graphics sites than this one, but I didn't want to bother looking for them. If all of the above is not self-evident, hit the reviews, Butikov (2004), "Square-wave excitation of a linear oscillator", American Journal of Physics 72  (4) 469-476.
