clear all;
close all;
s = 0.01; % an individual displacement of a molecule as a result of a single collision
n = input('Enter the number of molecular collisions, n = ');
t = 1 : n; % the vector of time instances
N = 500; 

A = rand(N,n) - 0.5; % A is an N-by-n matrix of random numbers equally distributed in the interval [-0.5, 0.5]
A = s*sign(A);
X = zeros(N,n); % X is an N-by-n matrix of total displacements of molecules 
                % columns of X correspond to positions of all N molecules at a time instance
for j = 1:n
        if ( j ~= 1)
            X(:,j) = sum(A(:,1:j)')';
            % the sum is taken for each row of A between 1 and j columns
            % the result is assigned to the j column of matrix X
            % X contains trajectories of all N molecules simultaneously
        else
            X(:,1) = A(:,1); 
            % no summations are performed for the first collision
        end    
end
    
Xav = sum(X)/N; % Xav is a row vector of n entries which give the mean at different time instances
XX = zeros(N,n);
for j = 1 : n
        XX(:,j) = Xav(j)*ones(N,1); % the j-th column of XX has the value of Xav(j)
end
Xav2 = sum((X - XX).^2)/N; % Xav2 is a row vector of n entries which give the variance at different time instances


figure(1);
plot(Xav); % plot the mean as function of time at figure(1)
title('Evolution of the mean \mu in time');
xlabel('Time'); ylabel('\mu'); 
figure(2);
plot(Xav2); % plot the variance as function of time at figure(2)
title('Evolution of the variance \sigma^2 in time');
xlabel('Time'); ylabel('\sigma^2'); 

% computations of the coefficient of diffusion D: sigma^2 = D*n