The Laplace transform is a useful tool for dealing with linear systems described by ODEs. As mentioned in another answer, the Laplace transform is defined for a larger class of functions than the related Fourier transform.
The 'big deal' is that the differential operator ('$\frac{d}{dt}$' or '$\frac{d}{dx}$') is converted into multiplication by '$s$', so differential equations become algebraic equations. In other words, convolution in the time or space-domain becomes multiplication in the s-domain. Another, often unspoken, 'big deal' is that the transform is unique in some sense (eg, if the transforms of two continuous functions agree, then the functions agree in the original domain). So if you can solve the problem in the s-domain, then you have solved it, in some sense, in the original domain. There is a formula for inversion, although tables are typically used for inversion. However, the inversion formula shows how the poles of the transformed functions manifest themselves in the time or space domain.
The Laplace transform comes in a few varieties; for engineering applications the most usual is the unilateral transform (behavior for $t<0$ is not relevant). Fourier transforms are often used to solve boundary value problems, Laplace transforms are often used to solve initial condition problems. Also, the Laplace transform succinctly captures input/output behavior or systems described by linear ODEs.
Regarding the 'frequency view'; instead of thinking of frequency as the $\omega$ in $\sin \omega t$, think of it as a collection of points in $\mathbb{C}$ that characterizes the behavior of $\hat{f} = \mathscr{L} f$.
For example, look at the Laplace transform of $f(t) = e^{\alpha t}$, which is $\hat{f}(s) = \frac{1}{s - \alpha}$. The single point $\alpha$ (which may be complex) completely characterizes the time domain behavior. More generally, the poles and zeros of $\hat{f}$ characterize the time domain behavior of $f$. Very loosely speaking, if $\hat{f}$ has poles $p_1,...,p_n$, then we expect $f$ to have time-domain 'behaviors' of the forms $e^{p_1 t},...,e^{p_n t}$ (the zeros, and pole multiplicities of $\hat{f}$ complicate this simplistic viewpoint somewhat). So, think of the frequencies (ie, poles & zeros) as characterizing the structure of $\hat{f}$.
In your question, I think you meant the system $\ddot {x} + \omega_n^2 x = f(t)$. The unilateral transform gives
$$s^2 \hat{x}(s) -s x(0) - x'(0)+ \omega_n^2 \hat{x}(s) = \hat{f}(s),$$
where $\hat{x}, \hat{f}$ are the Laplace transforms of $x,f$ respectively.
This equation is typically written in the following form, which shows the relationship between the input $\hat{f}$, the (time) initial conditions $x(0), x'(0)$, and the output $\hat{x}$:
$$\hat{x}(s) = \frac{s x(0) + x'(0)}{s^2 + \omega_n^2} + \frac{\hat{f}(s)}{s^2 + \omega_n^2}.$$
We can see that the $\frac{1}{s^2 + \omega_n^2}$ term 'contributes' two poles (at $s = \pm i \omega_n$) lying on the imaginary axis to $\hat{x}$. So, we expect (at least) behaviors involving $t \mapsto \sin \omega_n t$ and $t \mapsto \cos \omega_n t$.
If we take $f = 0$, you can see (meaning look up a table of transforms) that the initial conditions translate into a time function of the form $x(t) = x(0) \cos \omega_n t + \frac{x'(0)}{\omega_n} \sin \omega_n t$. So, in this particular problem, the initial conditions 'remain' forever.
If we take the system to be at rest initially (ie, take the initial conditions to be zero), then we need to know $\hat{f}$ in order to compute $\hat{x}$. If we take $f(t) = e^{i \omega t}$ (admittedly not real, but easier to compute), we have $\hat{f}(s) = \frac{1}{s-i\omega}$, which gives $\hat{x}(s) = \frac{1}{(s-i\omega)(s^2 + \omega_n^2)}$. If we take $w \neq w_n$, then using a partial fraction expansion we can write $\hat{x}(s) = \frac{1}{\omega_n^2 - \omega^2}(\frac{1}{s-i\omega} - \frac{s+i \omega}{s^2 + \omega_n^2})$, which gives $x(t) = \frac{1}{\omega_n^2 - \omega^2}(e^{i \omega t} - \cos \omega_n t - i \frac{\omega}{\omega_n} \sin \omega_n t)$. If $w = w_n$, then we obtain $\hat{x}(s) = \frac{i}{2 \omega_n}(\frac{1}{s^2+\omega_n^2} - \frac{1}{(s-i\omega_n)^2} ) $, which corresponds to $x(t) = \frac{1}{2 w_n^2} \sin \omega_n t - \frac{i}{2 \omega_n} t e^{i \omega_n t}$. Notice the response in this case is unbounded, even though the input is bounded. The '$t$' term arises because of the pole of multiplicity 2 at $s = i \omega_n$.
