Short answer: $\tau$ is the typical time it takes for heat (energy) to be transported over the distance $L$.
I'll try to elaborate a bit on your analogy to particle diffusion.
For particle diffusion in one dimension, you may think of the particle as jumping around on the x-axis. Some times it jumps to the right, and some times to the left. The end result is that it typically takes $\tau = L^2/\mathcal D$ to cover the distance $L$, when the diffusion constant is $\mathcal D$. The diffusion constant is a measure of how large the jumps are (in fact, how large the variance of the jumps is).
But for heat transport you may instead think of a chain of beads on the x-axis. Each bead is wiggling around its spot on the axis, and the more it wiggles, the higher the temperature at that position. Every now and then a wiggling bead will smack its neighbor, and exchange some energy with it. Some times the energy transfers from left to right, and some times the transfer is from right to left. The "thermal diffusivity" $\kappa$ is a measure of how often the beads collide, and how willing they are to exchange their energy with each other. The end result is that $\tau = L^2/\kappa$ is the typical time it takes for a "packet" of energy to travel over the distance $L$.
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.
Best Answer
Since the heat equation for 1D time-dependent conduction is $$k\frac{\partial^2 T}{\partial x^2}=\rho c\frac{\partial T}{\partial t},$$
the hand-wavy way to derive the characteristic time is to replace partial differentials with finite differences and note that the thermal diffusivity $D=k/\rho c$:
$$D\frac{\Delta T}{(\Delta x)^2}\approx\frac{\Delta T}{\Delta t}$$ $$\Delta t\approx\frac{(\Delta x)^2}{D}$$
But note that $L^2/D$ is not an exact time for anything to occur—it's just an estimate. In your example, the bar never reaches its steady state; it only asymptotically approaches it. (In practice, random fluctuations would soon hide any difference.) For a given geometry and material, the progress in reaching this steady state (e.g., the middle of the object is now halfway to its steady-state temperature) will scale with $L^2/D$.
You can often assume that the process of reaching equilibrium is well on its way after time $L^2/D$ has passed and that the process is nearly complete after several time constants have passed. However, I emphasize again that this is just an approximation, one that depends on the nature of the boundary conditions (e.g., constant temperature vs. constant flux), for example.