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](https://i.sstatic.net/9wlXp.jpg)) 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](https://i.sstatic.net/UUqtC.png)

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;