My question here is very computational. My problem is in mathematical physics, so I want to ask the community what kind of software they use to do the following computation if there is any? Let $$L_{n}=-\frac{n+1}{s} (u+v)n\frac{\partial}{\partial n} +\sum_{j=1}^{\infty} p_j (n+j)\frac{\partial}{\partial p_{n+j}}+\sum_{i+j=n}ij\frac{\partial^2}{\partial p_i \partial p_j} $$
For $n,m\geq 0$ it satisfies following equation $$[L_m,L_n]=(m-n)L_{m+n} $$ where $[]$ denote the commutation bracket in Weyl algebra. That is $[\frac{\partial}{\partial p_{j}}, p_j]=1$ other wise all other combination commute. Notice that there infinitely many $p_i$.
You can do this type of computation using the Maple Physics package; it supports commutator algebras, algebraic differential operators, tensorial non-commutative operators, etc. An example of a similar kind of calculation is presented in "Quantum Runge-Lenz Vector and the Hydrogen Atom, the hidden SO(4) symmetry", a step-by-step demonstration departing from basic principles using the kind of algebras and manipulation with differential operators you are asking about; at the end of it, there is a link to a PDF file with the steps visible. That post can give you an idea of how to formulate your problem.