I came across the System of equations
$\begin{equation}
\left\{ \begin{aligned}
\frac{d}{dt}p &= - \omega^2q \\
\frac{d}{dt}q &= p \\ p(0)&=p_0,q(0)=q_0
\end{aligned} \right.
\end{equation}$
which is equivalent to the standard form
$\begin{equation}
\left\{ \begin{aligned}
\frac{d}{dt}x(t) &= J^{-1}Ax(t) \\
x(0)=x0
\end{aligned} \right.
\end{equation}$ ,
where $$x=\begin{pmatrix}
p \\
q
\end{pmatrix} J^{-1}=\begin{pmatrix}
0 & -1 \\
1 & 0
\end{pmatrix}, A=\begin{pmatrix}
1 & 0 \\
0 & \omega^2
\end{pmatrix}$$
Now let $S=\begin{pmatrix}
\sqrt{\omega} & 0\\
0 & \dfrac{1}{\sqrt{\omega}}
\end{pmatrix}$, then $J^{-1}A$ gets transformed to $\tilde{A}=S^{-1}J^{-1}AS=\begin{pmatrix}
0 & -\omega \\
\omega & 0
\end{pmatrix}$
and the system of equations to
$\begin{equation}
\left\{ \begin{aligned}
\frac{d}{dt}\tilde{p} &= - \omega \tilde{q} \\
\frac{d}{dt}\tilde{q} &= \omega \tilde{p}
\\ \tilde{p}(0)&=\tilde{p}_0,\tilde{q}(0)=\tilde{q}_0
\end{aligned} \right.
\end{equation}$
The solution is given by $$\tilde{x}(t)=e^{t\tilde{A}}\tilde{x_0}=e^{S^{-1}J^{-1}AS}\tilde{x_0}=\begin{pmatrix} \cos(\omega t) & -\sin(\omega t)\\ \sin(\omega t) & \cos(\omega t) \end{pmatrix}\tilde{x_0}$$ Then I have to conclude $$x(t)= \begin{pmatrix} \cos(\omega t) & - \omega \sin(\omega t)\\ - \omega^{-1} \cos(\omega t) & \sin(\omega t) \end{pmatrix}x_0$$ which I think is an easy linear algebra exercise which unfortunately I can't solve and since I don't know which theory is involved, I can't rely on google either.. It is not clear to me how one computes $(\tilde{p},\tilde{q})$ in the first place, consequently I don't see how to invert this transformation. I tried by "doing the opposite", i.e. computing $S\tilde{x}S^{-1}$ but this doesn't yield the correct result. I'm not looking for the computations, hints on how to proceed are perfectly fine.
It seems your change of coordinates is $\tilde{x}(t) = S^{-1}x(t)$, because then we see $$ \frac{d\tilde{x}}{dt} = S^{-1}\frac{dx}{dt} = S^{-1}J^{-1}Ax = \tilde{A}S^{-1}x = \tilde{A}\tilde{x}. $$ You've computed $\tilde{x}(t)$ correctly from there. If we let $$ K = \begin{bmatrix} \cos(\omega t) & -\sin(\omega t)\\ \sin(\omega t) & \cos(\omega t)\end{bmatrix} $$ then we see $$ x(t) = S\tilde{x}(t) = SK\tilde{x}_0 = SKS^{-1}x_0. $$ So because we have $$ SKS^{-1} = \begin{bmatrix} \cos(\omega t) & -\omega\sin(\omega t)\\ \omega^{-1}\sin(\omega t) & \cos(\omega t)\end{bmatrix}, \text{ we have } x(t) = \begin{bmatrix} \cos(\omega t) & -\omega\sin(\omega t)\\ \omega^{-1}\sin(\omega t) & \cos(\omega t)\end{bmatrix}x_0 $$ or, in terms of coordinates, $$ p(t) = p_0\cos(\omega t) - \omega q_0\sin(\omega t),$$ $$q(t) = q_0\cos(\omega t) + \omega^{-1}p_0\sin(\omega t). $$ Hence there appears to be an error with your proposed solution. You can verify that this is the solution to your original system.
On your confusion about "how one computes $(\tilde{p}, \tilde{q})$ in the first place," I'm not sure what you mean. If you knew how to generate $S$ in the first place, then you ought to know that all you're doing is changing bases, hence the formula $\tilde{x} = S^{-1}x$. But if you mean that you don't know where $S$ came from, it arises from the attempt at "moving a copy of $\omega$" from the original $dp/dt$ equation to the $dq/dt$.
This is a kind of "symplectic scaling" which preserves the Hamiltonian structure of the equations. Basically, you need to scale $p$ by some factor $c$, but since there are only two coordinates being used, the only way to preserve the structure is to also scale $q$ by $1/c$. Play around with possible $c$ choices, and you'll arrive at your $S$ matrix.