p = 0.5;

%% (b)
nwalks = [100 1000, 10000];
for i=1:length(nwalks)
    subplot( 3, 1, i );
    for k=1:10
        x = randomwalk( nwalks(i), p );
        plot( x );
        hold on;
    end
    text( 0.05, 0.8, sprintf( 'N=%d', nwalks(i)), 'units', 'normalized' )
    xlabel( 'Number of steps' );
    ylabel( 'Position of walker' )
    hold off;
end
pause( 5.0 )

nsteps = 100;
nwalks = 10000;
subplot( 1, 1, 1 )
y = zeros( nwalks, nsteps );
for k = 1:nwalks
    x = randomwalk( nsteps, p );
    y(k,:) = x;    % store random walk
end
ns = 1:nsteps;
mu = mean(y, 1);
sdev = std(y, 1);
plot( ns, mu, 'b', 'linewidth', 4 )
hold on
plot( ns, sdev, 'r', 'linewidth', 4 )
xx = 0:0.01:nsteps;
plot( xx, sqrt(xx), 'k' )
plot( xx, zeros(length(xx),1), 'k' )
legend( 'mean', 'std', 'theory' )
hold off
pause( 3.0 );

%% (d) histograms:
tinx = [100, 30, 10];
colors = [ 0 0 1; 0.5 0 0.5; 1 0 0 ];
for i = 1:length(tinx)
    [h,b] = hist( y(:,tinx(i)), 20);
    h = h/sum(h)/(b(2)-b(1));
    bar(b, h, 1.0, 'facecolor', colors(i,:))
    hold on;
end
hold off;
xlabel('Position of walker');
ylabel('Probability density');