Stochastic dynamics of the QHO

This is the simplest way of simulating the stochastic evolution of a QHO using the Master equation discussed in class. 

We input:
 - The initial state (Subscript[n, 0])
 - The bath’s Bose-Einstein occupation (OverBar[N])
 - The time step Δt, which should be sufficiently small. 
 - The total number of time steps
 
 Output is written in a vector nv.
  
 This code is not meant to be efficient, but pedagogical. 
 The logic is as follows: if the system is in n, then it may either jump to n+1 or n-1 with probability 
 
(Subscript[P, up ]) = Δt (Subscript[W, n+1,n]) = Δt (OverBar[N]) (n+1)      and   (Subscript[P, down ]) = Δt (Subscript[W, n-1,n]) = Δt ((OverBar[N])+1)  n
 
 or it may stay where it is with probability (Subscript[P, stay]) = 1 - (Subscript[P, up])- (Subscript[P, down])
 This is simulated in Mathematica using the RandomChoice[] function.

Δt = 0.002;
Nb = 4.0;
n0 = 7;
steps = 6000;

Clear[nv];
nv = ConstantArray[0,steps]; (* Create a vector of empty entries *)
nv[[1]] = n0;
Do[
Pup = Δt Nb (nv[[i]]+1);
Pdown = Δt (Nb+1) nv[[i]];
Pstay = 1 - Pup - Pdown;
nv[[i+1]]=RandomChoice[{Pup,Pdown,Pstay}->{nv[[i]]+1,nv[[i]]-1,nv[[i]]}];

,{i,1,steps-1}];

ListLinePlot[nv,
(*GridLines->{None,{3}},*)
AspectRatio->1/5,
FrameTicks->{{Automatic,Automatic},{None,None}},
ImageSize->900,
BaseStyle->{FontFamily->"Times",20},
Frame->True,
Axes->None,
FrameLabel->{Style["time",Italic],Style["n",Italic]},
PlotRangePadding->None,
PlotStyle->red,
GridLines->{None,{Nb}},
GridLinesStyle->Directive[Gray, Dashed]
]