Great question. Of course, transforms are the best way to solve this, which you already know. But maybe when you start with a bunch of numbers from an oscilloscope instead of a nice, neat formula, your perspective changes a bit. Here is another way to approach the problem numerically. It is really best suited if you have some a priori info to bring to the problem, but that is often the case for such problems.
In the time domain, the output is \$y(t)=\int_{0}^tx(t')h(t-t')dt'\$. In this case, you know y(t) and x(t) and want to know h(t-t'). You can treat it as a very important class of problems called Inverse Problems. Now, if you knew nothing at all about \$h\$, this might be quite a bit more difficult, but if you can narrow down the space of solutions a bit, you have a good chance. The idea is to parameterize a potential impulse-response function, and then optimize the solution over those parameters. If you can come up with an initial guess to launch the optimization calculation, you have a good chance, but no guarantee, of success.
Since LTI systems so often have solutions such as exponentials and decaying oscillations, it gives you a basis for an educated guess. For instance, in this case, you see that the output y(t) involves an oscillatory component with frequency \$\omega = 4\$, so I don't think it requires any great leap to guess a second-order system solution of a decaying exponential with frequency \$\omega=4\$, such as \$h_o(t)=Ce^{-at}cos(\omega t + \phi)\$, with \$C\$, \$a\$,\$\phi\$, and maybe \$\omega\$ as optimization parameters.
Next you numerically calculate an output \$z(t)\$ to the known input \$x(t)\$, using a starting guess for the optimization parameters, and compute a mean-square error \$\int (y(t)-z(t))^2dt\$. You run a non-linear optimization over the parameters \$C\$, \$a\$, and \$\phi\$ (and maybe \$\omega\$) to get the solution that minimizes this error.
There is no guarantee you come up with a good solution this way, and it is dependent upon coming up with a decent 'guess' parameterization, but for many problems, this is not so unreasonable.
For more complicated system response functions, such as higher-order systems, your chances of solving the optimization problem diminish. But this is actually not all that difficult a method to implement, and could succeed for a good number of problems.
Of course, all this described as such is really nothing more than fitting data to a potential solution. However, looking at it as the more general inverse problem that it is motivates the use of other techniques from inverse problems to help deal with noise and overfitting, such as regularization.
