Problem: solve the EOM
$$\ddot x + \beta \dot x + \omega^2 x = 0$$$$\ddot x + \beta \dot x + \omega^2 x = f(t)$$
As an approach we shall use, additionally to $x(t), \dot x(t)$, two new parameters $y(t), \dot y(t)$.
Let us, magically, introduce a Lagrangian for this auxiliary system
$$L(x, y, \dot x, \dot y) = \dot x \dot y - \omega^2 x y - \beta \dot x y$$$$L(x, y, \dot x, \dot y, t) = \dot x \dot y - \beta \dot x y - \omega^2 x y - (x + y) f(t)$$
The important thing to notice is that the equations of motion for this system are
$$ \frac{d}{dt}\left( \frac{\partial L}{\partial \dot x} \right) - \frac{\partial L}{\partial x} = \ddot y - \beta \dot y + w^2 y = 0\\ \frac{d}{dt}\left( \frac{\partial L}{\partial \dot y} \right) - \frac{\partial L}{\partial y} = \ddot x + \beta \dot x + w^2 x = 0 $$$$ \frac{d}{dt}\left( \frac{\partial L}{\partial \dot x} \right) - \frac{\partial L}{\partial x} = \ddot y - \beta \dot y + w^2 y - f(t) = 0\\ \frac{d}{dt}\left( \frac{\partial L}{\partial \dot y} \right) - \frac{\partial L}{\partial y} = \ddot x + \beta \dot x + w^2 x - f(t) = 0 $$
As one can see, we recover the equations of motion for our original system along with an auxiliary EOM.
From here on out, everything goes according to theory for Hamiltonian mechanics. We can find the generalized momenta:
$$ p_x = \frac{\partial L}{\partial \dot x} = \dot y - \beta y\\ p_y = \frac{\partial L}{\partial \dot y} = \dot x$$
And rewriting the Langrangian as a Hamiltonian
$$H(x, y, p_x, p_y) = p_x p_y + \omega^2 x y + \beta y p_y$$$$H(x, y, p_x, p_y, t) = p_x p_y + \omega^2 x y + \beta y p_y + (x + y) f(t)$$
The method is a bit more general, see Conservative perturbation theory for nonconservative systems which introduced me to the idea of auxiliary parameters by the example of the Van der Pol oscillator.
As far as I can see, this method should play nicely even when $x \in \mathbb R^n$ in which case you would also choose $y \in \mathbb R^n$.