Will the energy of a randomly driven harmonic oscillator grow to infinity or oscillate about a finite value?

If you have a Gaussian random force, the equation becomes a Langevin equation. In physicists notation, you would write $$ \ddot x + x = \lambda \xi(t) \tag{1}$$ with $\lambda$ the strength of the random force and (the bracket denotes the stochastic average) $$\langle \xi(t) \rangle = 0, \quad \text{and} \quad \langle \xi(t) \xi(t') \rangle = \delta(t-t')\,. \tag{2}$$ Note that in mathematics instead stochastic differential equations are more conventional.

Let us assume that $x(0)=\dot x(0)=0$. We can solve (1) to obtain $$ x(t) = \lambda \int_0^t\!\sin(t-t') \xi(t')\,dt'\,. \tag{3}$$ This is a stochastic solution as it depends on the random function $\xi(t)$. However, from (3) together with (2) we can calculate statistical predictions. For example the average position is given by $$\langle x(t) \rangle =0\,,$$ which is not unexpected (just compare it to a random walk). So on average the oscillator does diverge as it does not even move.

Of course, the more reasonable measure if the harmonic oscillator performs an unbounded oscillation is the variance. We obtain $$\langle x(t)^2 \rangle = \lambda^2 \int_0^t \int_0^t\!\sin(t-t') \sin(t-t'') \langle\xi(t')\xi(t'')\rangle\,dt''\,dt' =\lambda^2 \int_0^t \sin^2(t-t')\,dt' = \lambda^2 \left(\frac{t}2 - \frac{\sin(2t)}{4}\right)\,. $$

From this we see that the typical amplitude of the oscillation, given by $\sqrt{\langle x(t)^2 \rangle }$ behaves as $$ \sqrt{\langle x(t)^2 \rangle} \sim \lambda \sqrt{\frac{t}{2}}$$ for $t\to\infty$; i.e., the oscillation grows without bounds. However, the amplitude of the oscillation only grows as $\sqrt{t}$ instead of proportional to $t$.

Does anyone know of a simple function I could replace f with to generate a continuous white noise driver?

I'm not sure if this approach will help you, but here it is. My approach to problems such as the one you presented is related to implementation of said problem and then derive some, at least numeric, result. Having this said, I would do as follows.

Consider a finite discrete Gaussian white noise which is stored in an array $G(n\Delta t)$, with $n$ from $0$ to $N$. This can be interpolated using a polynomial $P(t)$ of order $N-1$ which is unique. This polynomial should be the solution you are looking for for $t\in[0,N\Delta t]$. From this you can compute your driver function using

$$f(t)=\frac{d^2 P}{dt^2}(t) + P(t)$$