Consider a set of N balls that start at the origin. In a given unit of time, $\delta t$, the balls have a probability $p = 0.5$ of jumping a distance $\delta x$ to the right, and the same probability of jumping $\delta x$ to the left.

It can be shown (see Ockendon et al) that the concentration, $c(x, t)$, that is, the number of balls per unit length, is given by solving the heat equation,

\begin{equation} \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} \end{equation}

where the diffusion coefficient, $D = \left(\frac{p \delta x^2}{\delta t}\right)$.

I wrote a very simple simulation in Matlab that does the stochastic procedure for $N = 1000$, $\delta t = 0.1$, $\delta x = 0.005$, and runs until reaching $t = 600$. After the simulation is done, the domain is split into bins centered at $0.015k$, $k = 0, \pm 1, \ldots$ and the balls tallied.

On top of this, the fundamental (Direc delta) solution of the heat equation is plotted. This has solution \begin{equation} u = \frac{A}{\sqrt{4 \pi D t}} \exp\left( - \frac{x^2}{4 D t}\right). \end{equation}

The problem is that I'm not sure what to use for $A$. I assumed that if $u$ is concentration per unit length, then I should have \begin{equation} \text{Number of balls in bin} = u(x, t) \cdot N \cdot \text{bin size}. \end{equation}

However, this seems to underpredict the results at the center. Does anybody know where I went wrong, or whether this is expected? I could not get better results by varying $N$. An image can be seen below:

Image link here; I'm unable to post images

The Matlab code is included if you want to run it yourself.

clear; close all

Nball = 1000;               % Number of balls
dt = 0.1; dx = 5e-3;        % time step, space step
Ntime = floor(60*10/dt);    % number of time steps
p = 0.5;                    % probability of left/right jump

u = zeros(Nball, 1);        % Start at zero

for j = 2:Ntime
    r = rand(1, Nball);
    
    tmp = r > (1-p/2);
    u(tmp) = u(tmp) + dx;   % Jump right
    
    tmp = r <= p/2;
    u(tmp) = u(tmp) - dx;   % Jump left
end

bin = 30;                   % divide domain into intervals bin*dx
rightdom = [0:bin*dx:2];
[n, xout] = hist(u(:,end), [-fliplr(rightdom(2:end)), rightdom]);

bar(xout, n, 'c');
hold on
    
x = linspace(-2, 2, 200);
A = Nball*(bin*dx);         % Is this the right amplitude?

D = p*dx^2/dt;
t = dt*Ntime;
v = A/sqrt(4*pi*D*t)*exp(-x.^2/(4*D*t)); % Gaussian
    
p = plot(x, v, 'r', 'LineWidth', 2);
xlabel('x'); ylabel('Number of balls per compartment');
legend('Stochastic', 'Heat Equation');
hold off;