Bootstrapping the (Subscript[𝒲, 3]) algebra

Setting up the stage by invoking the OPEdefs package and modifying a notation:

ClearAll;
Needs["OPEdefs`"]

<<Notation`

Notation[(NotationTemplateTag[𝟙]) ⟺ (NotationTemplateTag[One])]
SetOPEOptions[SeriesArguments,{z,0}]

The Virasoro algebra works as usual:

Bosonic[T];
OPE[T,T]=MakeOPE[{0,c/2 One, 0, 2T, T'}];

OPE[T,T]//OPEToSeries

Following Zamolodchikov’s steps, we introduce a new primary W of weight (Subscript[h, W])=3:

Bosonic[W];
OPE[T,W]=MakeOPE[{0,0, 0, 3 W, W'}];

OPE[T,W]//OPEToSeries

What about the OPE W × W? The maximum weight of operators appearing in the pole structure must be (Subscript[h, max])=5, and the lightest weight is the identity’s weight (Subscript[h, 𝟙])=0. 
({{(1/(z^6)), 𝟙}, {(1/(z^5)), -}, {(1/(z^4)), T}, {(1/(z^3)), W,∂T}, {(1/(z^2)), ∂W, ∂∂T, 𝒩(TT)=(Subscript[Λ, TT])}, {(1/z), ∂∂W, ∂∂∂T, ∂𝒩(TT)=(Subscript[Λ, TT])', 𝒩(TW)=(Subscript[Λ, TW])}})

A couple of quasi-primary normal orderings appear. These take the following shapes:
𝒩(TT) = N(TT)-(3/10)∂∂T
𝒩(TW) = N(TW)-(3/14)∂∂W

𝒩[T,T]=NO[T,T]-(3/10)T'';
𝒩[T,W]=NO[T,W]-(3/14)W'';

𝒩[T,T]'

And I will define some rules for the sake of good notation:

NormalOrder={NO[T,T]->(Subscript[Λ, TT])+(3/10)T'', NO[T,W]->(Subscript[Λ, TW])+(3/14)W'', NO[T',T]->(1/2)(Subscript[Λ, TT])'+(( Derivative[3][T])/15) };

OPE[W,W]=MakeOPE[{
(Subscript[a, 1]) One, 
0, 
(Subscript[a, 2]) T, 
(Subscript[a, 3]) W+ (Subscript[a, 4]) T', 
(Subscript[a, 5]) W'+(Subscript[a, 6]) T''+(Subscript[a, 7])𝒩[T,T], 
(Subscript[a, 8]) W''+(Subscript[a, 9]) T'''+(Subscript[a, 10]) 𝒩[T,W]+(Subscript[a, 11])𝒩[T,T]'
}];

OPE[W,W]/.NormalOrder//FullSimplify//OPEToSeries

Now we want to close the algebra. The problem of closing the algebra is converted through Jacobi relations into the problem of solving a system of equations:

List1=OPEJacobi[T,T,T]//Flatten//DeleteCases[#,0]&//GetCoefficients//Flatten;
List2=OPEJacobi[T,T,W]//Flatten//DeleteCases[#,0]&//GetCoefficients//Flatten;
List3=OPEJacobi[T,W,W]//Flatten//DeleteCases[#,0]&//GetCoefficients//Flatten;
List4=OPEJacobi[W,W,W]//Flatten//DeleteCases[#,0]&//GetCoefficients//Flatten;
System=Join[List1, List2,List3, List4];

Bootstrap=Solve[System==0, {(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[a, 4]),(Subscript[a, 5]),(Subscript[a, 6]),(Subscript[a, 7]),(Subscript[a, 8]),(Subscript[a, 9]),(Subscript[a, 10]),(Subscript[a, 11])}]

Finally we choose the same normalization of W as Zamolodchikov:

(Subscript[a, 1])=(c/3);

The final outcome of the bootstrap procedure is:

OPE[W,W]/.Bootstrap[[1]]/.NormalOrder//FullSimplify//OPEToSeries

This coincides with the original Zamolodchikov’s result: