How to solve this coupled linear ODE I've tried using mathematica and for some reason it's not working
I think it's a pretty standard problem but for some reason I'm having a hard time
This is essentially a variant of the Rabbi problem I'll really appreciate the help :)
I will appreciate both explanation on how to solve this using mathematica or a manual solution
\begin{align}
i\hbar \frac{d\alpha}{dt}&=\frac{\hbar\omega_0}{2}\alpha+\frac{dE}{2}e^{i\omega t}\beta\\
i\hbar \frac{d\beta}{dt}&=-\frac{\hbar\omega_0}{2}\beta+\frac{dE}{2}e^{-i\omega t}\alpha
\end{align}
 A: We can rewrite the system as
$$
\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}\alpha\\\beta\end{pmatrix}
=\begin{pmatrix}-i\omega_0/2&-\frac{idE}{2\hbar}e^{i\omega t}\\-\frac{idE}{2\hbar}e^{-i\omega t}&i\omega_0/2\end{pmatrix}
\begin{pmatrix}\alpha\\\beta\end{pmatrix}
$$
and let us introduce $\tilde\alpha=e^{i\omega_0 t/2}\alpha$ and $\tilde\beta=e^{-i\omega_0 t/2}\beta$, giving
$$
\frac{\mathrm{d}}{\mathrm{d}t}\begin{pmatrix}\tilde\alpha\\\tilde\beta\end{pmatrix}
=\underbrace{-\frac{idE}{2\hbar}\begin{pmatrix}&-e^{i(\omega+\omega_0) t}\\e^{-i(\omega+\omega_0) t}&\end{pmatrix}}_{A(t)}\begin{pmatrix}\tilde\alpha\\\tilde\beta\end{pmatrix}
$$
the $A$ here is much nicer to handle.
So we compute
\begin{align*}
[A]\text{-terms}&:\int_0^t A(\tau)\,\mathrm{d}\tau \\
&= \left(\frac{dE}{2\hbar(\omega+\omega_0)}\right)\begin{pmatrix}&e^{i(\omega+\omega_0)t}-1\\e^{-i(\omega+\omega_0)t}-1&\end{pmatrix}\\
[A^2]\text{-terms}&:\int_0^t A(\tau_1)\int_0^{\tau_1} A(\tau_2)\,\mathrm{d}\tau_2\,\mathrm{d}\tau_1\\
&= \left(\frac{dE}{2\hbar(\omega+\omega_0)}\right)^2
\begin{pmatrix}i(\omega+\omega_0)t+1-e^{i(\omega+\omega_0)t}\\&e^{-i(\omega+\omega_0) t}-1- i(\omega+\omega_0)t\end{pmatrix}\\
&\dots
\end{align*}
You should see that the $[A^n]$-term has two entries which are $\pm(dE/(2\hbar(\omega+\omega_0)))^n$ times the tail part of the series $e^{\pm i(\omega+\omega_0)t}$.  Hence summing gives ultimately
$$
I+\int_0^t A(\tau)\,\mathrm{d}\tau+\int_0^t A(\tau_1)\int_0^{\tau_1} A(\tau_2)\,\mathrm{d}\tau_2\,\mathrm{d}\tau_1+\dots = 
$$
$$
\begin{pmatrix}
\frac1{2\omega_1} ( \omega_+ e^{i\omega_- t/2} - \omega_- e^{i\omega_+ t/2})
&
\ast
\\
\ast
&
\frac1{2\omega_1}( \omega_+ e^{i\omega_1 t/2} - \omega_- e^{-i\omega_1 t/2} )
\end{pmatrix}
$$
where
$\omega_1=\sqrt{(\omega+\omega_0)^2 + 4 (\frac{dE}{2\hbar})^2}$, $\omega_\pm= \omega+\omega_0 \pm \omega_1$ and I left two entries ($\ast$) for you to calculate.
A: I think you are trying to solve simultaneous linear differential equations of first order of the form $$\frac{dx}{dt}=ax+by \\ \frac{dy}{dt}=cx+dy.$$ I would highly recommend to go through the book "An Introduction to Dynamical system and chaos" (Springer) by G. C. Layek.
A: I'll begin by recasting this problem in a more concretely physical way. Introducing $\Psi=(\alpha,\beta)^\top$, we have
$$i\hbar\frac{d\Psi}{dt}=i\hbar \frac{d}{dt}\begin{pmatrix} \alpha\\ \beta\end{pmatrix}=\begin{pmatrix}\hbar \omega_0/2 & e^{i\omega t}dE/2\\ e^{-i\omega t}dE/2 & -\hbar\omega_0/2\end{pmatrix}\begin{pmatrix} \alpha\\ \beta\end{pmatrix}=H(t)\Psi$$
which is the time-dependent Schrodinger equation with Hamiltonian
$$\displaystyle H(t)
=
\begin{pmatrix}\hbar \omega_0/2 & e^{i\omega t}dE/2\\ e^{-i\omega t}dE/2 & -\hbar\omega_0/2\end{pmatrix}.$$
The Hamiltonian can be expressed in terms of Pauli matrices as
\begin{align}
H(t)&=\frac{\hbar\omega_0}{2}\begin{pmatrix} 1& 0 \\ 0 & -1\end{pmatrix}+\frac{dE}{2}\cos\omega t\begin{pmatrix} 0 & 1 \\ 1 & 0\end{pmatrix}-\frac{dE}{2}\sin\omega t\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}\\
&=\frac{\hbar \omega_0}{2}\sigma_z+\frac{dE}{2}(\cos(\omega t)\sigma_x-\sin(\omega t)\sigma_y)
\end{align}
Note that the $dE$ portion amounts to rotating the operator $\sigma_x$ around the $z$-axis by an angle $-\omega$. Explicitly,
\begin{align}
e^{i\omega t\sigma_z/2}\sigma_x e^{-i\omega t \sigma_z/2}
&=\begin{pmatrix} e^{i\omega t/2} & 0 \\ 0 & e^{-i \omega t/2}\end{pmatrix} \begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix}
\begin{pmatrix} e^{-i\omega t/2} & 0 \\ 0 & e^{i\omega t/2}\end{pmatrix}\\
&=\begin{pmatrix} 0 & e^{i\omega t} \\ e^{-i\omega t} & 0\end{pmatrix}
\end{align}
and therefore the Hamiltonian may be expressed as
$$H(t)=\frac{\hbar\omega_0}{2}+\frac{dE}{2}e^{i\omega t\sigma_z/2}\sigma_x e^{-i\omega t \sigma_z/2}$$
This amounts to an electron subject to an external magnetic field $\vec{B}=\frac12 (dE\cos\omega t,dE\sin\omega t,\hbar\omega_0)$ which rotates in some cone around the $z$-axis.
This suggests choosing a reference frame which rotates along with the magnetic field. To do this, we introduce $\Psi_r =e^{-i\omega t\sigma_z/2}\Psi$ ($r$ for rotated). Then the relevant Schrodinger equation is
\begin{align}
i\hbar \frac{d}{dt}\Psi_r
&=i\hbar \frac{d}{dt}e^{-i\omega t\sigma_z/2}\Psi\\
&=\frac{\hbar\omega}{2} \sigma_z \Psi_r+e^{-i\omega t\sigma_z/2}i\hbar \frac{d}{dt}\Psi\\
&=\frac{\hbar\omega}{2} \sigma_z \Psi_r+e^{-i\omega t\sigma_z/2}H(t)\Psi\\
&=\frac{\hbar\omega}{2} \sigma_z \Psi_r+e^{-i\omega t\sigma_z/2}H(t)e^{i\omega t\sigma_z/2}\Psi_r\\
&=H_r\Psi_r
\end{align}
where we now have the time-independent Hamiltonian $\displaystyle H_r=\frac{1}{2}\hbar(\omega+\omega_0) \sigma_z+\frac{dE}{2}\sigma_x$
in the rotating frame. This Hamiltonian is solvable by hand. To facilitate this, let $$ \Omega=\frac12\sqrt{(dE/\hbar)^2+(\omega+\omega_0)^2},\quad \tan\phi = \frac{dE}{\hbar(\omega+\omega_0)},$$
with which we may write $H_r=\hbar\Omega [\cos(\phi)\sigma_z+\sin (\phi)\sigma_x].$ The same trick as before means that we may diagonalize the Hamiltonian as
$$H_r = e^{-i\phi \sigma_y/2}\hbar \Omega \sigma_z e^{i\phi\sigma_y/2}.$$
We now solve the Schrodinger equation as $\Psi_r(t)=U_r(t)\Psi_r(t=0)$ where
\begin{align}
U_r(t)=e^{-i H_r t/\hbar}
&=e^{-i\phi \sigma_y/2}e^{-i \Omega \sigma_z t}e^{i\phi\sigma_y/2}\\
&=\begin{pmatrix} \cos\frac{\phi}{2} & -\sin \frac{\phi}{2}\\ \sin \frac{\phi}{2} & \cos\frac{\phi}{2}\end{pmatrix}
\begin{pmatrix}e^{-i B_r/\hbar} & 0 \\ 0 & e^{i B_r/\hbar}\end{pmatrix}
\begin{pmatrix} \cos\frac{\phi}{2} & \sin \frac{\phi}{2}\\ -\sin \frac{\phi}{2} & \cos\frac{\phi}{2}\end{pmatrix}\\
&=\begin{pmatrix}
\cos \Omega t -i \cos\phi \sin\Omega t &
-i \sin\phi \sin \Omega t\\
-i \sin\phi \sin\Omega t &
 \cos\Omega t+i \cos\phi \sin\Omega t
\end{pmatrix}
\end{align}
Returning to the original reference frame, the solution takes the form
$$\Psi(t)=e^{i\omega t\sigma_z/2}\Psi_r(t)=e^{i\omega t\sigma_z/2}U_r(t)\Psi_r(0)=U(t)\Psi(0)$$
where the time-evolution operator is now
\begin{align} 
U(t)
&= e^{i\omega t\sigma_z/2}U_r(t)\\
&=
\begin{pmatrix} e^{i\omega t/2} & 0 \\ 0 & e^{-i\omega t/2}\end{pmatrix}
\begin{pmatrix}
\cos \Omega t -i \cos\phi \sin\Omega t &
-i \sin\phi \sin \Omega t\\
-i \sin\phi \sin\Omega t &
 \cos\Omega t+i \cos\phi \sin\Omega t
\end{pmatrix}\\
&=\begin{pmatrix}
e^{i\omega t/2}(\cos \Omega t -i \cos\phi \sin\Omega t) &
-i e^{i\omega t/2}\sin\phi \sin \Omega t\\
-i e^{-i\omega t/2}\sin\phi \sin\Omega t &
 e^{-i\omega t/2}(\cos\Omega t+i \cos\phi \sin\Omega t)
\end{pmatrix}\\
&=\begin{pmatrix}
e^{i\omega t/2}(\cos \Omega t -i \Omega^{-1}(\omega+\omega_0) \sin\Omega t) &
-i e^{i\omega t/2}\Omega^{-1}(dE/\hbar) \sin \Omega t\\
-i e^{-i\omega t/2}\Omega^{-1}(dE/\hbar) \sin\Omega t &
 e^{-i\omega t/2}(\cos\Omega t+i \Omega^{-1}(\omega+\omega_0) \sin\Omega t)
\end{pmatrix}
\end{align}
This represents the full closed-form solution. As a sanity check, if $dE=0$ then $\Omega=\frac12( \omega+\omega_0)$ and thus
$$U(t)\to
\begin{pmatrix}
e^{i\omega t/2}e^{-i(\omega+\omega_0)t/2} & 0\\ 0 &
 e^{-i\omega t/2}e^{-i(\omega-\omega_0)t/2}
\end{pmatrix}
=
\begin{pmatrix}
e^{-i\omega_0 t/2} & 0\\ 0 &
 e^{-i\omega_0 t/2}
\end{pmatrix}
=e^{-i\omega_0 t \sigma_z/2}
$$
which agrees with $H(t)= \dfrac{\hbar\omega_0}{2}\sigma_z$ in this limit.
