It sounds like you're trying to solve the Langevin Equation. This is a model of Brownian motion where the particle experiences stochastic kicks at discrete time intervals. Your force, in this case, is a random variable you draw from a distribution each time step (instead of being given by an explicit formula).
For the simplest case, the "kicks" are generated by thermal noise, and have a gaussian distribution.
From the Wikipedia article linked above, the equation of motion is:
$$
m \frac{d^2 \vec{x}}{dt^2} = - \lambda \frac{d\vec{x}}{dt} + \vec{f(t)}
$$
Here the $\lambda$ term is from viscous friction with the fluid, and $\vec{f}$ is the random force. For a particle of radius $a$ in a fluid of dynamic viscosity $\eta$, then Stokes' Law says $\lambda = 6 \pi \eta a$.
If time is continuous, the correlation function for the kicks is:
$$
\langle f_i(t) f_j(t') \rangle = 2 \lambda k_B T \delta_{ij}\delta(t-t')
$$
Where $T$ is the temperature and $f_i$ is the $i$'th component of the force. To generate a force $\vec{f}(t)$ obeying this correlation function you may want to investigate gaussian processes.
Finally, there IS a relation between these (microphysical) constants and the diffusion coefficient $D$. It is called the Stokes-Einstein relation and was the first example of the Fluctuation-Dissipation Theorem. For the Langevin equation above, the relation is:
$$
D = \frac{k_B T}{\lambda} \ .
$$
Note: the Langevin equation is NOT an ordinary differential equation. It is a Stochastic differential equation, so the usual (e.g. Runge-Kutta) numerical methods won't apply.
the collision sets the time scale over which the particle can travel freely on average. Based on this picture, semiclassical, you can use the relation you quoted to calculate the Diffusion constant.
Best Answer
The diffusion constant of small spheres in a three dimensional fluid was found by Einstein $$ D = \frac{T}{6\pi\eta a} $$ where $\eta$ is the viscosity of the fluid, and $a$ is the radius of the particles. This is based on the Stokes result for the drag on a sphere, which does not immediately generalize to 2-d (sometimes called Stokes paradox). The Stokes paradox was resolved by Oseen. Using his result, I get $$ D = \frac{T}{4\pi\eta}\log\left(\frac{4\eta}{a\rho}\sqrt{\frac{2M}{T}}\right) $$ where $\rho$ is the density of the fluid, and $M$ is the mass of the disks.
Derivation: This is based on two ingedients. The first is the Langevin equation $$ \frac{d}{dt}\vec{p}=-\eta_D\vec{p}+\vec{\xi} $$ where $\eta_D$ is the drag force and $\xi$ is a stochastic force. Averaging this equation over the noise, and using $$ M\big\langle \big(\frac{d\vec{x}}{dt}\big)^2\big\rangle = dT $$ where $M$ is the mass of the particles, and $d$ is the number of dimensions, I get $$ \langle \vec{x}^2\rangle = \frac{2dT}{m\eta_D}\, t. $$ Comparing to the mean square deviation from the diffusion equation $\partial_0n-D\nabla^2 n=0$ I get $$ D =\frac{T}{m\eta_D} $$ independent of $d$. This is the Einstein relation.
The second ingredient is to use the drag force predicted by Stokes drag. This is the drag on a sphere/disk in a viscous fluid in the limit of small Reynolds number. In 3d this is the well known result $$ \eta_D = \frac{6\pi\eta a}{M} $$ and in 2d it is the more subtle result (due to Oseen) $$ \eta_D = \frac{4\pi\eta}{M}\frac{1}{\log(4/Re)} $$ where $Re=ua\rho/\eta$. Here $u$ is the relative velocity of the disk and the fluid, and $\rho$ is the density of the fluid.
Comments: The result in 3d is simple and intuitive. $D\sim T/(\eta a)$, so diffusion is faster if the temperature is higher, the viscosity of the solvent is smaller, and the particles are smaller. The 2d result is unusual, because even though we find the same dependence on $T/\eta$, the diffusion constant is only weakly (logarithmically) dependent on the particle size $a$.
We could have guessed this on purely dimensional grounds. There is no formula for the drag force that scales as $\eta a$. Indeed, if instead of diffusion of macroscopic particles (experiencing drag) I study microscopic particles experiencing quantum mechanical scattering with scattering length $a_s$ I find a power law $D\sim \frac{T}{n}\frac{1}{\sqrt{mT}a_S^2}$ in 3d, and log dependence $D\sim \frac{T}{n}\log(\frac{1}{mTa_s^2})$ in 2d.
What is also peculiar is that the Langevin equation in 2d has logarithmic drag $p\log(p)$. I think that this is indeed ultimately related to a breakdown of the diffusion approximation in 2d. This is well known, see http://www.sciencedirect.com/science/article/pii/0378437178900419. This does not mean that the $\log(a)$ dependence is wrong, it means that there are long time tails in the correlation function.