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 = ');
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);
% A becomes a random matrix of individual displacements +s or -s: 
% rows of A correspond to N molecules 
% columns of A correspond to n collisions that happen with each of the N molecules simultaneously
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
plot(X'); % plot all individual trajectories simultaneously
title('Individual trajectories of molecules in Brownian motion');
xlabel('Time'); ylabel('Displacements of molecules'); hold off;
