(* starting at 0, (tau1 n^2) plays of game A, followed by (tau2 n^2) plays of game B *) n = 100; alpha = 1/4; tau1 = 12/5; tau2 = 12/5; For[lambda = 1, lambda <= 5, lambda++, Print[DateList[]]; Print[]; Print["lambda = ", lambda, ", ", "n = ", n]; rho = 1. - lambda/n; p[0] = rho^(1/alpha - 1)/(1 + rho^(1/alpha - 1)); p[1] = 1/(1 + rho); PB = ConstantArray[0, {2 (tau1 + tau2) n^2 + 1}]; (* This array keeps track of the distribution of the random walk after t steps, t<=(tau1+tau2)n^2. The possible values of the random variable are -(tau1+tau2)n^2,...,(tau1+tau2)n^2. Since 1- dim array entries are labeled from 1 to the dimension, we add (tau1+ tau2)n^2+1 to the value to get the array entry label. *) PB[[-1 + (tau1 + tau2) n^2 + 1]] = 1./2; PB[[1 + (tau1 + tau2) n^2 + 1]] = 1./2; (* the distribution after 1 step *) For[t = 2, t <= (tau1 + tau2) n^2, t++, For[i = -t, i <= t, i = i + 2, PB[[i + (tau1 + tau2) n^2 + 1]] = If [Abs[i - 1] <= t - 1, PB[[i - 1 + (tau1 + tau2) n^2 + 1]], 0] If[t <= tau1 n^2, 1./2, p[Boole[Mod[i - 1, 4 n] >= n]]] + If[Abs[i + 1] <= t - 1, PB[[i + 1 + (tau1 + tau2) n^2 + 1]], 0] If[t <= tau1 n^2, 1./2, 1 - p[Boole[Mod[i + 1, 4 n] >= n]]]]]; (* Notice that we are not overwriting the array in this recursive step. Instead entries in odd positions are used to compute entries in even positions or vice-versa. *) g[x_] := If[x < 1, x, (4 - x)/3]; Print[Show[ ListPlot[{Table[{2 i/n - (tau1 + tau2) n, (n/2) PB[[ 2 i + 1]]}, {i, (tau1 + tau2) n^2/2 - 3 n, (tau1 + tau2) n^2/2 + 3 n}], Table[{2 i/n - (tau1 + tau2) n, If[i >= (tau1 + tau2) n^2/2 - 3 n/2, (n/2) PB[[ 2 i + 1]], -0.1]}, {i, (tau1 + tau2) n^2/2 - 3 n, (tau1 + tau2) n^2/2 + 3 n}], Table[{2 i/n - (tau1 + tau2) n, If[i <= (tau1 + tau2) n^2/2 + n/2, (n/2) PB[[ 2 i + 1]], -0.1]}, {i, (tau1 + tau2) n^2/2 - 3 n, (tau1 + tau2) n^2/2 + 3 n}]}, Joined -> True, PlotRange -> {{-6.2, 6.2}, {0, 4}}, Filling -> {1 -> {{2}, Red}, 1 -> {{3}, Green}}, Axes -> {True, False}, PlotStyle -> {Black, Black, Black, Black}, GridLines -> {{-3, 1}, None}, GridLinesStyle -> Directive[Gray, Dashed, Thick], Ticks -> {Table[{k, If[Mod[k, 4] == 0 || Mod[k, 4] == 1, k, ""], {If[Mod[k, 4] == 0, .012, .006], 0}}, {k, -6, 6, 1}], Automatic}], Plot[g[Mod[x, 4]], {x, -6, 6}, PlotRange -> {0, 4}, Axes -> {True, False}, GridLines -> {{-3, 1}, None}, GridLinesStyle -> Directive[Gray, Dashed, Thick], PlotStyle -> Blue]]]; Print["sum of probabilities = ", Sum[PB[[2 i + 1]], {i, 0, (tau1 + tau2) n^2}]]; Print["areas of peaks = {", N[Sum[PB[[2 i + 1]], {i, 0, (tau1 + tau2) n^2/2 - 3 n/2 - 1}], 9], ",", N[Sum[ PB[[2 i + 1]], {i, (tau1 + tau2) n^2/2 - 3 n/2, (tau1 + tau2) n^2/2 + n/2 - 1}], 9], ",", N[Sum[PB[[ 2 i + 1]], {i, (tau1 + tau2) n^2/2 + n/2, (tau1 + tau2) n^2}], 9], "}"]; Print["heights of peaks = {", N[(n/2) PB[[(tau1 + tau2) n^2 - 4 n + 1]]], ",", N[(n/2) PB[[(tau1 + tau2) n^2 + 1]]], ",", N[(n/2) PB[[(tau1 + tau2) n^2 + 4 n + 1]]], "}"]; Print["mean displacement = ", Sum[PB[[1 + 2 i]] (-(tau1 + tau2) n^2 + 2 i)/n, {i, 0, (tau1 + tau2) n^2}]]];