The long-standing Gaussian product inequality (GPI) conjecture states that $E [\prod_{j=1}^{n}X_j^{2m_j}]\geq\prod_{j=1}^{n}E[X_j^{2m_j}]$ for any centered Gaussian random vector $(X_1,\dots,X_n)$ and $m_1,\dots,m_n\in\mathbb{N}$. In this paper, we describe a computational algorithm involving sums-of-squares representations of multivariate polynomials that can be used to resolve the GPI conjecture. To exhibit the power of this novel method, we apply it to prove two new GPIs: $E[X_1^{2m_1}X_2^{6}X_3^{4}]\ge E[X_1^{2m_1}]E[X_2^{6}]E[X_3^{4}]$ and $E[X_1^{2m_1}X_2^{2}X_3^{2}X_4^{2}]\ge E[X_1^{2m_1}]E[X_2^{2}]E[X_3^{2}]E[X_4^{2}]$.