na = 3 ;             (* number of oscillators in Einstein solid A *)
nb = 3 ;             (* number of oscillators in Einstein solid B *)
qtotal = 6 ;     (* total energy in system A + B *)

Ω1 = Table[       (* Calculate (Subscript[Ω, Total])((Subscript[q, A])) *)
Binomial[ na -1. + q , q ] Binomial[ nb - 1. + (6-q) ,6-q] ,
{ q,0,6 }
]

(* Plot (Subscript[Ω, Total])((Subscript[q, A])) vs (Subscript[q, A]) *)
ListPlot[ Ω1, Filling -> Axis, PlotRange -> {0,100}]

(* Now plot probability of observing macrostate with energy (Subscript[q, A]) in A *)totalMicrostates = Binomial[ na + nb - 1. + qtotal , qtotal ]

ListPlot[ Ω1/totalMicrostates , Filling -> Axis, PlotRange -> {0,0.25} ]

(* Now look at Einstein solid with bigger subsystems, more energy *)
na = 300  ;              (* Einstein solid A with na oscillators *)
nb = 200 ;              (* Einstein solid B with nb oscillators *)
qtotal = 100 ;      (* total amount of energy in A and B *)

Ω2 = Table[           (* multiplicity Ω((Subscript[q, A])) of energy q in A, qtotal-q in B *)
Binomial[ na + q - 1., q ]  Binomial[ nb + (qtotal-q)-1. ,qtotal-q ],
{ q,0,qtotal }
] ;

Take[ Ω2, 20 ]    (* display first 20 elements of list e2 *)

Print[ Max[ Ω2 ], "  ", Min[ Ω2 ], "  ", Max[ Ω2 ] / Min[ Ω2 ] ]

(* Plot (Subscript[Ω, Total])((Subscript[q, A])) vs (Subscript[q, A]) *)
ListPlot[ Ω2, Filling -> Axis, PlotRange -> All]

(* Now calculate probability of observing system A+B with energy (Subscript[q, A]) in A *)
totalMicrostates = Sum[ Ω2[[i]], {i,Length[Ω2]}] // N

probs = Ω2 / totalMicrostates  ; (* prob of observing energy q in A, qtotal-q in B *)
Take[ Ω2, 20 ]

Print[ Min[probs], "  ", Max[probs], "  ", Max[probs]/Min[probs] ]

ListLinePlot[ probs , Filling-> Axis, PlotRange -> All]

stirling[ n_ ] := Sqrt[2.π n] ((n/E))^n ;

Plot[ (stirling[x] - x!)/x!, {x,1,30}, PlotRange -> {0,-0.04}]