% diffusion on a grid by a random walk

% create a vector of positions
clear all
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

 
N = 1000; % number of objects; change this number to get more or less error
X = zeros(N,1); % initial position is zero
alpha = 0.1; %probability of movement to the left
beta = alpha;  %probability of movement to the right

M = 100;  % total number of time steps
m = [0:M];
XS = zeros(N,M+1);% stores the positions at each time step
YS = zeros(N,2*M); 
c = [alpha,alpha+beta,1]; %this is the reaction vector to determine the actual move

% Specify the movement vector
xm(1) = 1; %to the right
xm(2) = -1; %to the left
xm(3) = 0; % don't move 

for j = 1:M % for each of M steps
    R = rand(N,1); % generate a random number for each particle
    for k = 1:N
       ndx= min(find(R(k)<=c)); % determine how to move
       X(k) = X(k) + xm(ndx); %update the position for particle k
    end
       XS(:,j+1) = X; % store the answer
       
end
% create a time plot of the position with jumps
for j = 1:M
    YS(:,2*j-1) = XS(:,j);
    YS(:,2*j) = XS(:,j);
    T(2*j-1) = j-1;
    T(2*j) = j;
end
mx = max(max(YS));
mn = min(min(YS));

max = max(mx,abs(mn));
% plot kk = 20 of these trajectories
kk = 20;
figure(1)
plot(T,YS(1:kk,:)','linewidth',2)
xlabel('time step','fontsize',20)
ylabel('position','fontsize',20)
axis([0 T(end) -max max])
  
 
figure(2)
%create a histogram of final positions

 Xtot = max-mn + 1;
Xlast = zeros(Xtot,1);
for j = 1:N
    pos = X(j) -mn+1;
    Xlast(pos) = Xlast(pos) + 1;
end
p = [mn:1:max];
plot(p,Xlast,'*')

% and compare this with the normal distribution
% Caution:  This is correct only if alpha = beta
gaussian = exp(-p.^2/(4*alpha*M))/sqrt(4*pi*M*alpha);

plot(p,Xlast/N,'*',p,gaussian,'--','linewidth',2)
xlabel('Position','fontsize',20)
ylabel('Probability','fontsize',20)

% plot the mean squared displacement
% Caution: if alpha ne beta, this is not the variance
V=sum(XS.^2)/N;
 
%
al_est = sum(m.*V)/sum(m.^2) % estimate for the slope of the variance

figure(3)
plot(m,V,'*',m,2*alpha*m,'--','linewidth',2)
ylabel('Mean Squared Displacement','fontsize',20)
xlabel('time step','fontsize',20)
  
