tic
L = 0.025;
tau = 0.0010;
NT = 100;             % Maximum steps
N = 1000;              % Number of particles to be considered
% variables related to generate many epsilons
dim3 = 30;                 %How many different epsilons to consider.
halfdim3 = dim3 / 2;
epsamp = 0.4;
epsstep = epsamp / dim3;       

for j = 0 : dim3
    % current epsilon, goes from -0.2 .. 0.2
    % Positive means tendency to go right, negative to left
    % for eps = 0.05, p = 0.55 and q = 0.45
    eps = epsstep * (j - halfdim3);
    % 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 a matrix, where each row
    % represents a particle and each column represents a timestep.
    R = cumsum(sign(-0.5 + eps + 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.

    % Let's transform the original matrix a bit - subtract X
    for i = 1 : N
        R(i, :) = R(i, :) - X;
    end
    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
    RESULT1(j + 1, :) = X;              %results of <X>
    RESULT2(j + 1, :) = X2;             %results of <(X-<X>)^2>
end
toc
% Let's plot all results in 3D this time..
[XX, YY] = meshgrid(tau:tau:NT*tau, -halfdim3 * epsstep:epsstep:halfdim3*epsstep);
figure(1)
surf(XX, YY, RESULT1);
xlabel('Time in s')
ylabel('epsilon')
zlabel('value of X')
title('X dependence of time and epsilon. Dependence of epsilon really exists in this graph. Dependence of time is linear.');
view([5,5])
colorbar
figure(2)
surf(XX, YY, RESULT2);
xlabel('Time in s')
ylabel('epsilon')
zlabel('value of <(x - <x>)^2>')
title(' <(x - <x>)^2> dependence of time and epsilon. Dependence of epsilon does not really exist in this graph. Dependence of time is linear');
view([-45,50])
colorbar