This system can be solved exactly by Floquet theory. The Floquet theory deals with periodically driven quantum systems, and this is exactly the prototypical example where a solution can be obtained.
Floquet theory in a nutshell
Consider the time-dependent Schrodinger equation with periodic time-dependent Hamiltonian $H(t)$ (such that $H(t+nT)=H(t)$, where $n$ is an integer and $T$ the period):
$$
i\hbar \partial_t \psi_{\alpha}(t) = H(t) \psi_{\alpha}(t). \;\;\;\;\;\;\;\;\; (1)
$$
In analogy with Bloch theory of particles in periodic potentials, we can assume that the eigenmode $\psi_{\alpha}(t)$ can be written as
$$
\psi_{\alpha}(t) = e^{-i\varepsilon_{\alpha}t/\hbar} u_{\alpha}(t) = \sum_{m=-\infty}^{+\infty} e^{-i(\varepsilon_{\alpha}/\hbar + m\omega)t} u_{\alpha}^{(m)}
\;\;\;\;\;\;\;\;\; (2)
$$
where $u_{\alpha}(t+nT) = u_{\alpha}(t)$ is a time periodic function with period T, and hence it can be expanded in a Fourier series with standard frequency $\omega = 2\pi/T$ and weights $u_{\alpha}^{(m)}$.
Substituting (2) into (1), you get the important equation
$$
\sum_n (H_{mn} - n\hbar\omega \delta_{mn}) u_{\alpha}^{(n)} = \varepsilon_{\alpha}u_{\alpha}^{(m)};
\;\;\;\;\;\;
H_{mn}\equiv \frac{1}{T}\int_0^T dt H(t) e^{i(m-n)\omega t}.
\;\;\;\;\;\;
(3)
$$
This can be recast into an eigenvalue problem by organizing the Fourier components into a vector: for example, in the two level system ($\alpha=1,2$), the vector $\vec{u}$ is:
$$
\vec{u} = ( \dots , u^{(-1)}_1, \,u^{(-1)}_2 , u^{(0)}_1, \,u^{(0)}_2 , u^{(1)}_1, \,u^{(1)}_2 , \dots)
$$
and the matrix we want to diagonalize is an infinite matrix made by $2\times 2$ blocks, each block being labeled by a pair of integers $(m,n)$.
Driven two level system
Let's now consider the following time periodic Hamiltonian that describes a two level system with a very specific periodic drive:
$$
H(t) = \frac{\varepsilon}{2} \, \hat{\sigma}_z + \frac{V}{2} \left( \cos(\omega t) \hat{\sigma}_x + \sin(\omega t) \hat{\sigma}_y \right),
$$
where $\hat{\sigma}_a$ is the $a$-th Pauli matrix. Notice that this is not the only possible choice for the drive, but this is the only one for which an analytic solution can be provided (as far as I know). The physical system we have in mind here is a spin-$1/2$ in a magnetic field along $z$ and a circularly polarized magnetic field (with constant modulus) in the $xy$ plane.
The matrix to be diagonalized here is the following:
\begin{equation}
\frac{1}{2} \left(
\begin{array}{c|cc|cc|cc|c}
\ddots & & & & & & & \\
\hline
& \varepsilon+2\hbar\omega & 0 & 0 & 0 & 0 & 0 & \\
& 0 & -\varepsilon+2\hbar\omega & V & 0 & 0 & 0 & \\
\hline
& 0 & V & \varepsilon & 0 & 0 & 0 & \\
& 0 & 0 & 0 & -\varepsilon & V & 0 & \\
\hline
& 0 & 0 & 0 & V & \varepsilon - 2\hbar\omega & 0 & \\
& 0 & 0 & 0 & 0 & 0 & -\varepsilon - 2\hbar\omega & \\
\hline
& & & & & & & \ddots
\end{array}
\right),
\end{equation}
where each box is a $2\times 2$ block labeled by $(m,n)$. You can see immediately that this matrix is block diagonal (just remove the lines and look :) ), so you can easily find the eigenvalues:
$$
\varepsilon_{\pm, n} = \frac{\hbar \omega \pm \hbar \omega_R}{2} - n\omega,
\;\;\;\;\;\;\;\;\;
(\hbar\omega_R)^2 \equiv (\hbar\omega - \varepsilon)^2 + V^2
$$
and the corresponding eigenvectors ($\left|\uparrow\right\rangle = (1,0)^T \;\; \left|\downarrow\right\rangle = (0,1)^T$):
$$
u_{\pm}^m = \frac{ V \delta_{m,0} \left|\downarrow\right\rangle + (\varepsilon - \hbar\omega \pm \hbar \omega_R) \delta_{m,+1} \left|\uparrow\right\rangle }{\sqrt{V^2 + (\varepsilon - \hbar\omega \pm \hbar \omega_R)^2}}.
$$
From these solutions you can get important dynamical properties of the system, e.g. the transition probability to go from $\left|\downarrow\right\rangle$ at $t=0$ to $\left|\uparrow\right\rangle$ at time $t$. I will be very sketchy here: the idea is to write down the wave function as a superposition of the two eigenmodes $\psi(t) = a \psi_+(t) + b \psi_-(t)$, find the coefficients $a$ and $b$ by the condition $\left\langle\downarrow \right| \psi(0) \rangle = 1$ and finally computing
$$
P(t) \equiv |\left\langle \uparrow \right| \psi(t) \rangle|^2 = \frac{V^2}{\omega_R^2} \sin^2{(\frac{\omega_R t}{2})}
$$
