function randwalk(p)

N = 500;
J = 25;
L = 1;

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

figure;
x_sq = zeros(J,N);
x_avg = zeros(J,1);
x_f = zeros(N,1);
for n=1:N
  fprintf('n=%d/%d\n',n,N);

  x = zeros(J,1);
  for j=2:J
    k = 2*(rand > p)-1;
    x(j) = x(j-1) + k*L;

    subplot(2,2,1)
    plot([0:j-1],x(1:j));
    title(sprintf('random walk path\n(%d/%d steps taken)',j,J));
    axis([0 J-1 -2*sqrt(J) 2*sqrt(J)])
    xlabel('number of steps');
    ylabel('x');

    drawnow

  end
  hold on;
  plot(J-1,x(end),'rx');
  hold off;

  x_f(n) = x(end);

  x_avg = (x_avg*(n-1) + x) ./ n;

  subplot(2,2,2);
  plot([0:J-1],x_avg);
  hold on;
  plot([0 J-1],[0 0],'r--')
  hold off;
  title(sprintf('average distance\n(%d/%d paths in average)',n,N));
  axis([0 J-1 -2*sqrt(J) 2*sqrt(J)])
  xlabel('number of steps');
  ylabel('<x>')

  x_sq(:,n) = x.^2;

  x_rms = sqrt(mean(x_sq(:,1:n),2));

  subplot(2,2,3)
  hold off;
  plot([0:J-1],x_rms);
  hold on;
  plot([0:J-1],sqrt([0:J-1]),'r--')
  title(sprintf('rms distance\n(%d/%d paths in average)',n,N));
  axis([0 J-1 0 2*sqrt(J)])
  xlabel('number of steps');
  ylabel('{<x^2>}^{1/2}')

  subplot(2,2,4)
  hist(x_f(1:n))
  title(sprintf('distribution of final positions\nafter %d steps',J));
  axis(round([-3*sqrt(J) 3*sqrt(J) 0 N/3]))
  xlabel('x_{25}')
  ylabel('count')
  h = findobj(gca,'Type','patch');
  set(h,'FaceColor','r')

  drawnow

end

