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 molecules in an ensemble, N = ');
n = input('Enter the number of molecular collisions, n = ');
t = 1 : n; % the vector of time instances

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
Dtime = Xav2./t; % the value of the coefficient diffusion D depends on the time t, when we compute it
D = sum(Dtime)/n; % this is a naive and rough way to approximate the theoretical value of D
                     % NOTE: a correct solution of the approximation problem should be done with 
                     
plot(t,X','b'); % plot all individual trajectories simultaneously
hold on;
x1 = sqrt(D*t); % the theoretical time dependence of the standard deviation displays 
                % a typical trajectory of an "average" molecule
x2 = -x1; % due to the symmetry, an "average" molecule can travel both to the right and to the left
plot(t,x1,'r',t,x2,'r');
title('Trajectories of molecules of an ensemble and the trajectories of an average molecule');
xlabel('Time'); ylabel('Displacements of molecules'); hold off;
