function diffuse1d(N,T,D)

if ~exist('N')
  N = 1e4;
end

if ~exist('T')
  T = 1e3;
end

if ~exist('D')
  D = 0.5;
end

xmax = 4*sqrt(2*D*T);
nmax = N/2;

x = zeros(N,1);
xrms = zeros(N,1);

figure;
for t=1:T
  xrms(t) = std(x);

  subplot(1,3,1)
  plot(x,zeros(size(x)),'rx');
  axis([-xmax xmax -1 1]);
  xlabel('position')
  ylabel('(meaningless)');
  title(sprintf('particle position\n(%d particles)',N));

  subplot(1,3,2);
  hist(x);
  axis([-xmax xmax 0 nmax])
  title(sprintf('distribution of particle positions\n(%d of %d steps)',t,T));
  xlabel('position');
  ylabel('count');

  subplot(1,3,3);
  plot(xrms(1:t),'x');
  hold on;
  plot(sqrt(2*D*[1:T]),'r--');
  hold off;
  axis([0 T 0 sqrt(1.25*2*D*T)]);
  title('standard deviation of particle position over time');
  xlabel('time');
  ylabel('standard deviation');
  legend('actual','predicted');


  drawnow;


  x = x + 2*(rand(size(x)) > .5) - 1;
end