The exact solution of your equation can be written as

$$x(t) = x(0) + \frac{m}{6 \pi a \nu} \dot{x}(0) - \frac{m}{6 \pi a \nu} \dot{x}(0) \, \text{e}^{-\frac{6\pi a\nu}{m}t} + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} F(\tau_2) \, . $$

$x(0)$ and $\dot{x}(0)$ are the initial conditions. Then whatever you want to compute, you can write it as an expression that depends on $F(t)$ and take its average over the fluctuations of $F(t)$.

For example

\begin{align*} \langle x(t) \rangle = & \langle x(0) \rangle + \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle- \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle \, \text{e}^{-\frac{6\pi a\nu}{m}t} \\ & + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \langle F(\tau_2) \rangle \, .\end{align*}

If your particle starts at rest and at the origin,

$$ \langle x(0) \rangle = 0 \, \qquad \langle \dot{x}(0) \rangle = 0 \, ,$$

and if it experiences a constant force,

$$ \langle F(t) \rangle = F \, ,$$

then you find

\begin{align*}\langle x(t) \rangle & = \frac{F}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \\ & = \frac{F}{m} \frac{m}{6 \pi a \mu} \left[ t + \frac{m}{6 \pi a \nu} \left( \text{e}^{-\frac{6 \pi a \nu}{m}t}-1\right) \right] \, .\end{align*}

You see that you can choose the statistics of $F(t)$ freely. Then if you know the moments of $F(t)$, you can compute the moments of $x(t)$. It's all about computing integrals. Typically on chooses Gaussian statistics with

$$ \langle F(t) \rangle = 0 \, , \qquad \langle F(t_1) F(t_2) \rangle = D \, \delta(t_1-t_2) \, .$$

If you insist on solving this problem numerically, you need to discretise time

$$ t \in \left[0,\infty\right[ \rightarrow t \in \left\{t_i\right\}_{i = 1,..,N} \, . $$

Then you can sample $F(t)$ according to your favourite probability distribution

$$ P\left[F(t_1),F(t_2),..,F(t_N)\right] \, .$$

For every sample you get a discretised function, $F(t_i)$ and you can switch to finite difference derivatives (for example) to solve your differential equation. Then you average at the end.

The exact solution of your equation can be written as

$$x(t) = x(0) + \frac{m}{6 \pi a \nu} \dot{x}(0) - \frac{m}{6 \pi a \nu} \dot{x}(0) \, \text{e}^{-\frac{6\pi a\nu}{m}t} + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} F(\tau_2) \, . $$

$x(0)$ and $\dot{x}(0)$ are the initial conditions. Then whatever you want to compute, you can write it as an expression that depends on $F(t)$ and take its average over the fluctuations of $F(t)$.

For example

\begin{align*} \langle x(t) \rangle = & \langle x(0) \rangle + \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle- \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle \, \text{e}^{-\frac{6\pi a\nu}{m}t} \\ & + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \langle F(\tau_2) \rangle \, .\end{align*}

If your particle starts at rest and at the origin,

$$ \langle x(0) \rangle = 0 \, \qquad \langle \dot{x}(0) \rangle = 0 \, ,$$

and if it experiences a constant force,

$$ \langle F(t) \rangle = F \, ,$$

then you find

\begin{align*}\langle x(t) \rangle & = \frac{F}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \\ & = \frac{F}{m} \frac{m}{6 \pi a \mu} \left[ t + \frac{m}{6 \pi a \nu} \left( \text{e}^{-\frac{6 \pi a \nu}{m}t}-1\right) \right] \, .\end{align*}

You see that you can choose the statistics of $F(t)$ freely. Then if you know the moments of $F(t)$, you can compute the moments of $x(t)$. It's all about computing integrals. Typically one chooses Gaussian statistics with

$$ \langle F(t) \rangle = 0 \, , \qquad \langle F(t_1) F(t_2) \rangle = D \, \delta(t_1-t_2) \, .$$

If you insist on solving this problem numerically, you need to discretise time

$$ t \in \left[0,\infty\right[ \rightarrow t \in \left\{t_i\right\}_{i = 1,..,N} \, . $$

Then you can sample $F(t)$ according to your favourite probability distribution

$$ P\left[F(t_1),F(t_2),..,F(t_N)\right] \, .$$

For every sample you get a discretised function, $F(t_i)$ and you can switch to finite difference derivatives (for example) to solve your differential equation. Then you average at the end.

The exact solution of your equation can be written as

$$x(t) = x(0) + \frac{m}{6 \pi a \nu} \dot{x}(0) - \frac{m}{6 \pi a \nu} \dot{x}(0) \, \text{e}^{-\frac{6\pi a\nu}{m}t} + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} F(\tau_2) \, . $$

$x(0)$ and $\dot{x}(0)$ are the initial conditions. Then whatever you want to compute, you can write it as an expression that depends on $F(t)$ and take its average over the fluctuations of $F(t)$.

For example

\begin{align*} \langle x(t) \rangle = & \langle x(0) \rangle + \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle- \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle \, \text{e}^{-\frac{6\pi a\nu}{m}t} \\ & + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \langle F(\tau_2) \rangle \, .\end{align*}

If your particle starts at rest and at the origin,

$$ \langle x(0) \rangle = 0 \, \qquad \langle \dot{x}(0) \rangle = 0 \, ,$$

and if it experiences a constant force,

$$ \langle F(t) \rangle = F \, ,$$

then you find

\begin{align*}\langle x(t) \rangle & = \frac{F}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \\ & = \frac{F}{m} \frac{m}{6 \pi a \mu} \left[ t + \frac{m}{6 \pi a \nu} \left( \text{e}^{-\frac{6 \pi a \nu}{m}t}-1\right) \right] \, .\end{align*}

You see that you can choose the statistics of $F(t)$ freely. Then it's all about computing integrals. Typically on chooses Gaussian statistics with

$$ \langle F(t) \rangle = 0 \, , \qquad \langle F(t_1) F(t_2) \rangle = D \, \delta(t_1-t_2) \, .$$

The exact solution of your equation can be written as

$$x(t) = x(0) + \frac{m}{6 \pi a \nu} \dot{x}(0) - \frac{m}{6 \pi a \nu} \dot{x}(0) \, \text{e}^{-\frac{6\pi a\nu}{m}t} + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} F(\tau_2) \, . $$

$x(0)$ and $\dot{x}(0)$ are the initial conditions. Then whatever you want to compute, you can write it as an expression that depends on $F(t)$ and take its average over the fluctuations of $F(t)$.

For example

\begin{align*} \langle x(t) \rangle = & \langle x(0) \rangle + \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle- \frac{m}{6 \pi a \nu} \langle \dot{x}(0) \rangle \, \text{e}^{-\frac{6\pi a\nu}{m}t} \\ & + \frac{1}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \langle F(\tau_2) \rangle \, .\end{align*}

If your particle starts at rest and at the origin,

$$ \langle x(0) \rangle = 0 \, \qquad \langle \dot{x}(0) \rangle = 0 \, ,$$

and if it experiences a constant force,

$$ \langle F(t) \rangle = F \, ,$$

then you find

\begin{align*}\langle x(t) \rangle & = \frac{F}{m}\int_0^t \text{d}\tau_1 \int_0^{\tau_1} \text{d}\tau_2 \text{e}^{\frac{6\pi a\nu}{m} \left(\tau_2-\tau_1\right)} \\ & = \frac{F}{m} \frac{m}{6 \pi a \mu} \left[ t + \frac{m}{6 \pi a \nu} \left( \text{e}^{-\frac{6 \pi a \nu}{m}t}-1\right) \right] \, .\end{align*}

You see that you can choose the statistics of $F(t)$ freely. Then if you know the moments of $F(t)$, you can compute the moments of $x(t)$. It's all about computing integrals. Typically on chooses Gaussian statistics with

$$ \langle F(t) \rangle = 0 \, , \qquad \langle F(t_1) F(t_2) \rangle = D \, \delta(t_1-t_2) \, .$$

If you insist on solving this problem numerically, you need to discretise time

$$ t \in \left[0,\infty\right[ \rightarrow t \in \left\{t_i\right\}_{i = 1,..,N} \, . $$

Then you can sample $F(t)$ according to your favourite probability distribution

$$ P\left[F(t_1),F(t_2),..,F(t_N)\right] \, .$$

For every sample you get a discretised function, $F(t_i)$ and you can switch to finite difference derivatives (for example) to solve your differential equation. Then you average at the end.