How to code oscillator driven by Gaussian white noise? Edit: How to convert ODE to a system of SDE's?

What you are dealing with is called stochastic differential equation. Go back to the differential form: $$ \mathbf{X}_t = \left[\begin{array}{c} X_t \\ \dot{X}_t \end{array} \right],$$ and write down the equation in matrix form $$d\mathbf{X}_t = \mathbf{M} \cdot \mathbf{X}_t dt + \left[\begin{array}{c}dW_t\\0\end{array}\right],$$ where $dW_t = \eta(t)dt$ and $\mathbf{M} = \left[\begin{array}{cc}0 & 1 \\ -\omega_0^2 & 0\end{array} \right]$. Now you can numerically simulate the process using Euler-Maruyama method: $$\mathbf{X}_{t+1} = \mathbf{X}_t + \mathbf{M} \cdot \mathbf{X}_t \Delta t + \left[\begin{array}{c}\Delta W_t\\0\end{array}\right],$$ and keep in mind that $\Delta W_t$ is a Gaussian random variable (with parameters mentioned in the question). If your discretization domain is small enough and you collected enough samples you should see a plot similar to the one below. Blue line is the mean $\langle X_t\rangle$ and orange $\sqrt{\langle X_t^2 \rangle}$.

EDIT

A little bit of theoretical explanation. The solution can be written as $$\mathbf{X}_t = e^{t \mathbf{M}} \mathbf{X}_0 + \int\limits_{0}^{t} e^{-(s-t)\mathbf{M}} \left[\begin{array}{c}\eta(s)\\0\end{array}\right]ds.$$ Since $\mathbf{X}_0 = \mathbf{0}$, we can write $$X_t = \frac{1}{\omega_0}\int\limits_{0}^{t} \cos[\omega_0(s-t)] \eta(s)ds$$ and you get (for $\omega_0 = 1$) $$\langle X_t^2\rangle = \frac{1}{2}t + \frac{\sin(2t)}{4}$$

SAMPLE Python code

# -*- coding: utf-8 -*-
import numpy as np

def run(x0=np.array([.0,.0]), n=40000, dt=1.0E-04, omega=1.0):

    sol = np.array([])

    M = np.array([[0, 1.],[-omega**2, 0.]])

    x = x0.copy()
    for i in range(0,n):
        sol = np.append(sol, x[0])
        x += M @ x * dt + np.array([1.,0.]) * np.random.normal(scale=np.sqrt(dt))


    return sol

sol = np.array([run() for i in range(0,500)])

mean  = np.mean(sol, axis=0)
sigma = np.sqrt(np.var(sol, axis=0) + mean**2)

dt = 1.0E-04
x = np.arange(0, len(mean))
y = np.sqrt(x * dt/2. + np.sin(2. * x * dt)/4.)

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(15,10))


ax.plot(x, mean)
ax.plot(x, sigma)
ax.plot(x, y)
plt.show()