Inspired by some intrinsic relations between Coulomb gas integrals and Gaussian multiplicative chaos, this article introduces a general mechanism to prove BPZ equations of order $(r,1)$ and $(1,r)$ in the setting of probabilistic Liouville conformal field theory, a family of conformal field theory which depends on a parameter $γ\in (0,2)$. The method consists in regrouping singularities on the degenerate insertion, and transforming the proof into an algebraic problem. With this method we show that BPZ equations hold on the sphere for the parameter $γ\in [\sqrt{2},2)$ in the case $(r,1)$ and for $γ\in (0,2)$ in the case $(1,r)$. The same technique applies to the boundary Liouville field theory when the bulk cosmological constant $μ_{\text{bulk}} = 0$, where we prove BPZ equations of order $(r,1)$ and $(1,r)$ for $γ\in (0,2)$.