% diffusion on a grid by a random walk
%Here we allow alpha and beta to be different
% 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 right
beta = 0.2;  %probability of movement to the left

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;  %this can be a problem if pos exceeds the range of X definition
    Xlast(pos) = Xlast(pos) + 1;
end
p = [mn:1:max];
plot(p,Xlast,'*')

% and compare this with the normal distribution
% there is a drift term and the diffusion coefficient is modified
gaussian = exp(-(p-(alpha-beta)*M).^2/(2*(alpha+beta)*M))/sqrt(2*pi*M*(alpha+beta));

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

% find the mean of the displacement
mn = sum(XS)/N;
figure(3)
plot(m,mn,'*',m,(alpha-beta)*m,'--')
ylabel('Mean','fontsize',20)
xlabel('time step','fontsize',20)
 
%calculate and plot the variance
V=sum((XS-ones(N,1)*mn).^2)/N;
 
%
al_est = sum(m.*V)/sum(m.^2) % estimate for the slope of the variance

figure(4)
plot(m,V,'*',m,(alpha+beta)*m,'--','linewidth',2)
ylabel('Variance','fontsize',20)
xlabel('time step','fontsize',20)

