In physics, the simple harmonic oscillator with a time-dependent frequency $\omega(t)$ obeys the differential equation
\begin{align} \frac{d^2x(t)}{dt^2}+\omega(t)^2x(t)=0\ . \end{align}
If $\omega$ is a constant, then
\begin{align} x(t)=C_1e^{i\omega t}+C_2e^{-i\omega t}\ . \end{align}
How would I find a solution to the differential equation if $\omega(t)$ is an arbitrary time-dependent function? I tried looking in books on solutions to differential equations (such as Handbook of exact solutions for ordinary differential equations). All solutions I found are for a specific form of $\omega(t)$, rather than an arbitrary $\omega(t)$.
I also tried the substitution $x(t)=s(t)e^{i\gamma(t)}$. Assuming $s(t)$ and $\gamma(t)$ are both real functions, then subbing $x(t)$ into the above differential equation and splitting the result into its real and imaginary terms, I find the $s(t)$ and $\gamma(t)$ satisfy
\begin{align} \frac{d^2s(t)}{dt^2}-s(t)\left(\frac{d\gamma(t)}{dt}\right)^2+\omega(t)^2s(t)=0\ ,\\ s(t)\frac{d^2\gamma(t)}{dt^2}+2\frac{ds(t)}{dt}\frac{d\gamma(t)}{dt}=0\ , \end{align}
but these equations are more complicated than my starting point.
There is no general closed-form solution. For small $\omega(t)$, one can often obtain a power series approximation (or use the Method of Frobenius if a regular singular point occurs).
For large $\omega(t)$, an approximate solution can be generated using the Liouville–Green formulae: put $x(t) = e^{S(t)} $, and the equation becomes $$ S'^2 + S'' = \omega^2. $$ If you now assume a solution of the form $S = S_0 + S_1 + S_2 + \dotsb$, where $S_k \gg S_{k+1}$ this equation is formally equivalent to the equations $$ S'^2 = -\omega^2 \\ 2S_0' S_1' + S_0'' = 0 \\ 2S_{k}S_0 + S_{k-1}'' + \sum_{n=1}^{k-1} S_n' S_{n-k}' = 0, $$ which can be solved order-by-order. The first two terms are $S_0 = \pm i\int^t \omega(s) \,ds$ and $S_1 = -\frac{1}{2}\log{\omega}$.
There are conditions that are required for this to work, namely the above inequalities between the $S_k$, and $S_{N+1} \ll 1$ so that the error in the exponential is small ($e^{O(1)}$ is not $1+O(1)$ in general).
In particular, in many cases this gives $e^{S_0+S_1} = \omega^{-1/2}\exp{\left(\pm i\int^t \omega(s) \, ds \right)}$ as the leading-order behaviour in $\omega$, which does agree with the solution you have when $\omega$ is constant.