Alternative ways to perform the same task

I have a code that I used to simulate the diffusion of a particle. I have done it using for loops. My next step will involve creating four simulations (4 particles) and adding interaction terms. I believe that writing a code to do that using for loops will be get messy pretty fast. My plan to make things easier to do involves creating an array, where the cells will be matrices for each simulation. I can then use cellfun to manipulate the entries of each cell. However, I am concerned with a beginning step. Looking at the first cell, how would I go about calculating the mean (sum column / #columns), of the first cell. Below is my code, showing how I did this with a for loop.

mu = 0;
kb=physconst('Boltzmann');
sig = sqrt((2*kb*300));
dt = 1e-3;
t = 0:dt:1;
P=2000;% Time vector
x = zeros(length(t),P); % Allocate output vector, set initial conditionrandn
rng('Default')
rng(1);  

for p=1:P;
for i = 1:length(t)-1
    x(i+1,p) = x(i,p)+sig*sqrt(dt)*randn;
end
end
figure;
plot(t,x);
X=x.';
A=mean(X);
B=mean(X.^2);
figure;
plot(t,A);
figure;
plot(t,B)
slope=polyfit(t,B,1)

Answer 1

[Erik Keever]

I think I might first ask, "do I need to write new code to do this?" Both particle and particle-in-cell simulation codes for diffusion, DLA, plasma physics, Boltzmann-type equations, Smolchuski-type equations, Navier-Stokes and many more have been written, debugged, tested, optimized, parallelized and made available...

If the answer is still yes, first consider how you'd rewrite this in vector notation:

for p=1:P;
for i = 1:length(t)-1
    x(i+1,p) = x(i,p)+sig*sqrt(dt)*randn;
end
end

Or you'll likely find the execution time grinding to a halt pretty quickly (profile on/profile off/profview are friends). As to the particular question you're asking: "Looking at the first cell, how would I go about calculating the mean (sum column / #columns), of the first cell."... is there a reason

mean(theCellArray{n})

wouldn't work? You might also consider extending the array into the 3rd dimension rather than using cells, if equal numbers of particles of each species is acceptable. Most functions which don't explicitly depend on specific array shapes (i.e. other than linear algebra itself) will extend intuitively with 3- or higher-dimensional arrays. Consider also if it's necessary to store the -entire- time history of the entire simulation, since that has a way of eating up memory, or if a finite amount of total history can be stored and much smaller running outputs computed based on that.