Our aim is to write and execute a Monte Carlo program to model diffusion in one-dimension.
Assume that a system of N independent particles all begin at the position x = 0.
At each timestep, a particle can move with equal probability either
left or right a single step of 
.
The following is a Mathematica program which calculates and outputs the distribution of particles
after M steps.
Needs["Statistics`DataManipulation`"]
Needs["Statistics`DescriptiveStatistics`"]
ranstep := 2*Random[Integer,{0,1}] - 1
ranwalk[n_] := FoldList[Plus, 0, Table[ranstep, {n}]]
mrwlast[m_Integer,n_Integer] := Table[Last[ranwalk[n]], {i,m}]
pmrw[m_,n_,xmin_,xmax_,dx_] := BinCounts[mrwlast[m,n], {xmin, xmax, dx}]
ranwalk[100]
{0, 1, 0, -1, 0, -1, 0, -1, 0, -1, 0, -1, -2, -1, 0, -1, 0, -1, -2,
-1, 0, -1, 0, -1, -2, -3, -4, -5, -6, -5, -4, -3, -2, -1, -2, -3,
-4, -5, -4, -5, -6, -7, -6, -5, -4, -3, -4, -5, -4, -5, -6, -5,
-6, -7, -8, -9, -10, -9, -8, -7, -8, -9, -10, -11, -10, -11, -10,
-9, -10, -11, -12, -11, -10, -9, -10, -9, -8, -9, -8, -7, -8, -7,
-6, -7, -6, -7, -8, -9, -8, -9, -10, -9, -8, -9, -8, -7, -6, -7,
-6, -7, -6}
mrwall[m_Integer,n_Integer] := Table[ranwalk[n], {i,m}]
N[Table[StandardDeviation[Transpose[mrwall[100,100]] [[i]] ], {i,1,100}]]
{0, 1., 1.40576, 1.75821, 1.93845, 2.24508, 2.36336, 2.50446, 2.86702,
2.94351, 3.30112, 3.42987, 3.66314, 3.79132, 3.90369, 4.21996,
4.25054, 4.43471, 4.4687, 4.68292, 4.87484, 5.00743, 5.17867,
5.14748, 5.18366, 5.32344, 5.55974, 5.57542, 5.73784, 5.79442,
5.89127, 5.93224, 5.93214, 6.00996, 6.03354, 6.1065, 6.20912,
6.20798, 6.301, 6.2969, 6.52984, 6.53766, 6.45406, 6.56849, 6.70402,
6.9205, 6.91218, 7.00361, 6.87211, 7.07718, 6.9931, 7.09745, 7.04878,
7.07281, 7.08676, 7.14103, 7.30584, 7.21981, 7.11021, 7.1043,
7.16363, 7.21533, 7.33708, 7.32779, 7.58244, 7.54834, 7.58744,
7.73425, 7.73341, 7.7062, 7.69793, 7.56972, 7.60924, 7.42573,
7.51306, 7.46169, 7.64513, 7.75837, 7.70871, 7.99606, 7.84731,
7.84702, 7.88757, 7.94933, 7.98673, 8.02468, 8.16635, 8.26518,
8.50134, 8.58307, 8.81438, 9.06709, 9.11207, 9.16954, 9.24993,
9.24746, 9.369, 9.48524, 9.41907, 9.53674}
ListPlot[%]
-Graphics-
Fit[%55, {Sqrt[i]}, i]
0.940735 Sqrt[i]
ListPlot[Table[%58, {i,100}]]
-Graphics-
Show[%56,%60]
-Graphics-