tic
L = 0.025;
tau = 0.0010;
D = L^2 / tau
NT = 100;    % Maximum steps
N = 10000;  % Number of particles to be considered

%we create an N times NT matrix of random values, and then cumsum over
% all elements in the rows. So as a result, we have matrix, where each row
% represents a particle and each column represents a timestep.
R = cumsum(sign(-0.5 + rand(N, NT)), 2) * L;

T = (1 : NT);                   %array of steps, later array of times...
X(1 : NT) = sum(R(:,T)) / N;        %average over all rows? vertically I mean.
X2(1 : NT) = sum(R(:,T).^2) / N;    %average of squared values over all rows.
T = T * tau;                    %let's make T correct for plotting

%fit
fit = polyfit(T, X2, 1);     %slope is the first.

% as Slope = C*D, then C = Slope / D
C = fit(1) / D
toc

figure('Name', 'X dependence on time')
plot(T, X);
figure('Name', 'X^2 dependence on time, fitted line also...')
plot(T, X2, T, T*fit(1)+fit(2));
