%----------------------------------------


%   Tony Hyun Kim


%   18.354, PS#1 (Problem 3)


%----------------------------------------


function randomwalk





N = 10000; % Num particles


T = 10;   % Num steps





x = zeros(1,N);


lim  = floor(20);


bins = [-lim:.5:lim];





% Display properties


figure;


axis([-lim lim 0 1/3*N]);


xlabel( 'Position' );


ylabel( 'Num particles' );


title( ['Random walk: ' num2str(N) ' particles'] );


grid on;





for t = 1:T


    cla();


    hold on; 


    tic


    hist(x,bins); % Draw the random walkers


    rt = toc;


    drawsoln(t,rt,N,lim);


    hold off;


    drawnow;


    for n = 1:N


       x(n) = update(x(n)); 


    end


end


return





%----------------------------------------


% Walking rules


%----------------------------------------


function x_n1 = update(x_n)





delta = .5;


k = 2;     % Num outcomes 


d = mod( floor(k*rand()), k );


switch d


    case 0


        x_n1 = x_n + delta;


    case {1, 2, 3}


        x_n1 = x_n - delta;


end


return





%----------------------------------------


% Draw continuum soln.


%----------------------------------------


function drawsoln(t,rt, N, lim)





t = t*rt;


D = 0.5*(.5^2)/rt; 


x = linspace(-lim, lim, 300);


y = N/sqrt(4*D*pi*t)*exp(-x.^2/(4*D*t));


plot(x, y);





return