Is there another way to solve this differential equation?
One, not-entirely-standard solution (though it does generalise, to some extent):
Write your equation as $\frac{\ddot x}{x}=-\omega^2$, and let $u=\frac{\dot x}{x}$.
Then $\dot u=\frac{\ddot x}{x}-\left(\frac{\dot x}{x}\right)^2\implies \dot u+u^2+\omega^2=0\implies\frac{\omega\dot u}{u^2+\omega^2}+\omega=0$.
From this we see:
$$\dot{\left(\arctan\frac{u}{\omega}\right)}+\omega=0\implies u=-\omega\tan(\omega t+c)$$
for some $c$. Thus
$$\frac{\dot x}{x}=-\omega\tan(\omega t+c)\implies \dot{(\log x)}=\dot{(\log \cos (\omega t+c))}$$
and we deduce that $x=R\cos(\omega t+c)$ for some $R,c$.
Since your second order linear ordinary differential equation has constant coefficients, you can generalise the solution to your equation (below) for any values of $\gamma$ and $\omega$ without any limitations by using the concept of the characteristic polynomial.
$$\frac{d^2 z}{dt^2}+2\gamma\frac{dz}{dt}+\omega^2 z=0$$
We will assume the form of the solution to be $z(t)=e^{\lambda t}$
We will substitute this into our differential equation to obtain:
$$\lambda^2 e^{\lambda t}+2\gamma\lambda e^{\lambda t}+\omega^2 e^{\lambda t}=0$$
$$e^{\lambda t}(\lambda^2+2\gamma\lambda+\omega^2)=0$$
We deduce that $e^{\lambda t}$ will only equal zero for non-zero values of $t$ if $\lambda \to -\infty$. Hence, this is a trivial solution.
Thus, we must only consider the solution to $\lambda$ on $\lambda^2+2\gamma\lambda+\omega^2=0$ (Our characteristic polynomial).
We solve this using the quadratic formula. Therefore, the two roots of the characteristic polynomial will be:
$$\lambda=\frac{-2\gamma \pm \sqrt{4\gamma^2-4\omega^2}}{2}=-\gamma \pm \sqrt{\gamma^2-\omega^2}$$
We denote $\lambda_1$ and $\lambda_2$ to be the roots of the above equation.
We notice that depending on the values of the constants $\gamma$ and $\omega$, the solutions to $\lambda$ will be different, and so hence our general solutions will be significantly different.
We separate each case:
- If both roots are real and distinct (i.e. $\lambda_1 \neq \lambda_2$)
- If both roots are complex.
- Real repeated roots (i.e. $\lambda_1=\lambda_2$)
Case 1
Since the particlular solutions are $z_1(t)=e^{\lambda_1 t}$ and $z_2(t)=e^{\lambda_2 t}$, they will yield the general solution when substituted to $z(t)=k_1 z_1(t) + k_2 z_2(t)$.
$$z(t)=k_1 e^{\lambda_1 t}+k_2 e^{\lambda_2 t}$$
Case 2
We denote a complex number by $\lambda=a+bi$.
Since $\lambda_{1,2}$ will both be complex, using $e^{i \theta}=\cos{\theta}+i\sin{\theta}$, we obtain the general solution:
$$z(t)=k_1 e^{at} \cos{(bt)} + k_2 e^{at} \sin{(bt)}$$
Case 3
$$z(t)=k_1 e^{\lambda t}+k_2 t e^{\lambda t}$$
END CASES
Hence, for the solution to the differential equation you wanted to find, $\ddot{z}+\omega z=0$, let $\gamma=0$. Since both solutions to $\lambda$ are purely imaginary (i.e. Complex with $a=0$), the form of the solution will be Case 2. The general solution is hence:
$$z(t)=k_1 \cos{(\omega t)} + k_2 \sin{(\omega t)}$$
You can substitute your initial conditions into the equation to obtain the values of the arbitrary constants $k_1$ and $k_2$ for simple harmonic motion as:
$$z(t)=z_0 \cos{(\omega t)} + \frac{v_0}{\omega} \sin{(\omega t)}$$
Which can be written in the form:
$$z(t)=A (\cos{\omega t+\varphi})$$
For more information, see:
Real, Distinct Roots, Complex Roots and Repeated Roots
Let's first make clear that we are not assuming that the solution for this differential equation has exponential form. This equation belongs to a class of differential equations whose have simple solutions and indeed have exponential form on a more subtle sense. Now we will find out what class of differential equations are we talking about by defining a new variable.
Let $ v = \dot z $. Then $ \dot v = \ddot z = - 2\gamma\dot z - \omega_0^2z $.
So we can see that this equation is actually is a system of homogenous linear differential equations. Namely $$ \begin{cases} \dot z = v \\ \dot v = - 2\gamma v - \omega_0^2z \end{cases} $$
Which can be written in matrix form by setting $ x = \begin{bmatrix} z \\ v \end{bmatrix} $ and $ A = \begin{bmatrix} 0 & 1 \\ - \omega_0^2 & - 2\gamma \end{bmatrix} $.
Using this notation the system simply becomes $$ \dot x = Ax $$
This equation doesn't remember you anything? You should be familiar with this equation when $ A \in \mathbb{R} $ and its solution $ x(t) = ke^{At} $ where $k$ is some real constant. It is from here that arises the motivation for matrix exponentiation. Consider the following question:
It is possible to extend the definition exponentiation from real numbers to real square matrices such that it preserves all its properties including the function $ x(t) = ke^{At} $ being the solution of the (now matrix) equation $ \dot x = Ax $ ?
The answer is yes and it is done naturally by just subtituting real numbers by real matrices in the power series definition of the exponential function. That is
$$ e^A := \sum _{ k=0 }^{ \infty }{ \frac { 1 }{ k! } { A }^{ k } } $$
Even more astonishing is that this definition even works for complex matrices!
Now back to the particular case $ \ddot x +\omega_0^2 x=0 $, the reason its real solution is written in terms of trigonometric functions is that its correspoding matrix has complex eigenvalues and Euler's formula $$ e^{iz} = cos(z) + i sin(z) $$ plays a role here.
By using the following two propositions which can be proven true and Euler's Formula we derive the general real solution from the general complex solution.
Proposition 1. The general complex solution of the linear system $ \dot x = Ax $ where $A$ is a $n$ x $n$ real matrix with $n$ distinct (possible complex) eigenvalues $\lambda_1, ..., \lambda_n$ associated with eigenvectors $v_1, ..., v_n$ is given by $$ x(t) = \sum _{ i=1 }^{ n }{ c_i e^{\lambda_i t} v_i } $$ where $c_i$ are complex constants. Note that this is complex exponentiation.
Proposition 2. The general real solution of the linear system $ \dot x = Ax $ where $A$ is a $n$ x $n$ real matrix with $n$ distinct (possible complex) eigenvalues $\lambda_1, ..., \lambda_n$ associated with eigenvectors $v_1, ..., v_n$ is a linear combination of real and imaginary parts of complex solutions. That is, if $x(t)$ is a complex solution, then $y(t) = k_1 Re(x(t)) + k_2 Im(x(t))$ is a real solution, where $k_1$ and $k_2$ are real constants.
(The details are rough but this is supposed to be just a glimpse of the theory behind)
Solution. Writing the equation associated matrix $ A = \begin{bmatrix} 0 & 1 \\ - \omega_0^2 & 0 \end{bmatrix} $ we find its eigenvalues $ \lambda_1 = -i w_0 $ and $ \lambda_2 = i w_0 $ with corresponding eigenvectors $v_1$ and $v_2$. Now we apply (Prop. 1) to get the complex general solution, followed by Euler's formula to get complex linear combinations of trigonometric functions and finnaly (Prop. 2) with some algebra to arrive at the general real solution given by: $$ x(t) = \begin{bmatrix} z(t) \\ v(t) \end{bmatrix} = \begin{bmatrix} k_1 cos(\omega_0 t) + k_2 sin(\omega_0 t) \\ - k_1 sin(\omega_0 t) + k_2 cos(\omega_0 t) \end{bmatrix} $$
You have to go through a bit of theory but after it is done it gets much more easier to find the general solution for this class of differential equations.
I hope I brought up some insight of a natural way to dealing with this kind of equations. You can find the details, worked examples and more about in the great book entitled Differential Equations, Dynamical Systems, and an Introduction to Chaos by Hirsch, Smale and Devaney.