Yes, such formulas exist (they are related to the so-called "linearization formulas"), and can be found in the literature. For instance, for the standard Hermite polynomials $H_k(x)$, the corresponding formula is:
$$
\int_{-\infty}^\infty H_\ell(x) \, H_m(x) \, H_n(x) \, e^{-x^2} dx = \frac{2^{(\ell+m+n)/2}\ell! m! n! \sqrt{\pi}}{(\frac{\ell+m-n}{2})!(\frac{m+n-\ell}{2})!(\frac{n+\ell-m}{2})!}
$$
if $\ell+m+n$ is even and the sum of any two of $\ell,m,n$ is not less than the third, or zero otherwise. This is equation (6.8.3) in the book "Special Functions", by Andrews-Askey-Roy.
Notice that your weight function is slightly different (it's a scaled version, corresponding to the "probabilists' Hermite polynomials"), but an easy change of variables will give the result that you want.
Here is one approach to prove the orthogonality relation for the Legendre polynomials $P_n (x)$.
We wish to prove
\begin{align*}
\int^1_{-1} P_m (x) P_n (x) \, dx = \begin{cases}
0, & m \neq n\\[1ex]
\displaystyle{\frac{2}{2n + 1}}, & m = n
\end{cases}
\end{align*}
The case $m \neq n$
As $P_m (x)$ and $P_n (x)$ satisfy the Legendre differential equation, namely
$$(1 - x^2) \frac{d^2 P_m}{dx^2} - 2x \frac{dP_m}{dx} + m(m + 1) P_m =
0,$$
and
$$(1 - x^2) \frac{d^2 P_n}{dx^2} - 2x \frac{dP_n}{dx} + n(n + 1) P_n =
0,$$
each of these equations may be rewritten as
\begin{align}
\frac{d}{dx} \left [(1 - x^2) \frac{dP_m}{dx} \right ] + m(m + 1)
P_m &= 0, \qquad (1)
\end{align}
and
\begin{equation}
\frac{d}{dx} \left [(1 - x^2) \frac{dP_n}{dx} \right ] + n(n + 1)
P_n = 0. \qquad (2)
\end{equation}
Multiplying Eq. (1) by $P_n (x)$ and
Eq. (2) by $P_m(x)$ and subtracting we have
$$P_n \frac{d}{dx} \left [(1 - x^2) \frac{dP_m}{dx} \right ] -
P_m \frac{d}{dx} \left [(1 - x^2) \frac{dP_n}{dx} \right ] + [m(m +
1) - n(n + 1)]P_m P_n = 0. \,\, (*)$$
Now recognising
$$\frac{d}{dx} \left [P_n (1 - x^2) \frac{dP_m}{dx} \right ]
= \frac{dP_n}{dx} (1 - x^2) \frac{dP_m}{dx} + P_n \frac{d}{dx} \left
[(1 - x^2) \frac{dP_m}{dx} \right ],$$
Eq. ($*$) may be rewritten as
$$\frac{d}{dx} \left [P_n(1 - x^2) \frac{dP_m}{dx} \right ]
- \frac{dP_n}{dx} (1 - x^2) \frac{dP_m}{dx} - \frac{d}{dx} \left
[P_m(1 - x^2) \frac{dP_n}{dx} \right ] + \frac{dP_m}{dx} (1 -
x^2) \frac{dP_n}{dx} + (m^2 + m - n^2 - n) P_m P_n = 0,$$
which reduces to
$$\frac{d}{dx} \left [P_n(1 - x^2) \frac{dP_m}{dx} \right ]
- \frac{d}{dx} \left [P_m(1 - x^2) \frac{dP_n}{dx} \right ] + (m^2 + m -
n^2 - n) P_m P_n = 0.$$
On integrating up the above equation with respect to $x$ from $-1$ to 1 we have
\begin{align*}
\int^1_{-1} \left \{\frac{d}{dx} \left [P_n(1 -
x^2) \frac{dP_m}{dx} \right ]
- \frac{d}{dx} \left [P_m(1 - x^2) \frac{dP_n}{dx} \right ] \right \} dx\\
\hspace{-8.0cm}+ (m^2 + m - n^2 - n) \int^1_{-1} P_m P_n \, dx &= 0\\
\left [P_n(1 - x^2) \frac{dP_m}{dx} - P_m(1 -
x^2) \frac{dP_n}{dx} \right ]^1_{-1} + (m - n)(m + n +
1) \int^1_{-1} P_m P_n \, dx &= 0\\
(m - n)(m + n + 1) \int^1_{-1} P_m P_n \, dx &= 0.
\end{align*}
Note the term appearing in the square brackets vanishes at both the upper and lower limits due to the presence of the factor $(1 - x^2)$ .
So when $m \neq n$, the factor out the front of the integral cancels
and we obtain
$$\int^1_{-1} P_m (x) P_n (x) \, dx = 0,$$
as required to prove.
The case $m = n$
To prove the $m = n$ case let us define
$$A_n = \int^1_{-1} [P_n (x)]^2 \, dx.$$
From Bonnet's recurrence relation for the Legendre polynomials, namely
$$(n + 1) P_{n + 1} (x) - (2n + 1)x P_n (x) + n P_{n - 1} (x) = 0,$$
if we replace $n$ with $n - 1$, after rearranging we have
$$P_n (x) = \frac{2n - 1}{n} x P_{n - 1} (x) - \frac{n - 1}{n} P_{n -
2} (x).$$
Using this result, substituting for one of the $P_n(x)$ terms in the
integral for $A_n$ we have
\begin{align*}
A_n &= \int^1_{-1} P_n (x) P_n (x) \, dx\\
&= \int^1_{-1} P_n (x) \left [\frac{2n - 1}{n} x P_{n - 1} (x) - \frac{n - 1}{n} P_{n -
2} (x) \right ] \, dx\\
&= \frac{2n - 1}{n} \int^1_{-1} x P_n (x) P_{n - 1} (x) \, dx
- \frac{n - 1}{n} \int^1_{-1} P_n (x) P_{n - 2} (x) \, dx
\end{align*}
But by the orthogonality property, as
$$\int^1_{-1} P_n (x) P_{n - 2} (x) \, dx = 0,$$
this yields
$$A_n (x) = \frac{2n - 1}{n} \int^1_{-1} x P_n (x) P_{n - 1} (x) \, dx.$$
Now from Bonnet's recurrence relation, by making the term $x P_n (x)$ the
subject we have
$$x P_n (x) = \frac{1}{2n + 1} \left [(n + 1) P_{n + 1} (x) + n P_{n -
1} (x) \right ],$$
which on substituting into the integral for $A_n$ one obtains
\begin{align*}
A_n &= \frac{2n - 1}{n} \cdot \frac{n + 1}{2n +
1} \int^1_{-1} \!\! P_{n +
1} (x) P_{n -1} (x) dx + \frac{2n - 1}{2n + 1} \int^1_{-1} [P_{n -
1} (x)]^2 dx\\
&= \frac{2n - 1}{2n + 1} A_{n - 1}, \quad n = 1,2,3,\ldots
\end{align*}
where we have again made use of the orthogonality property.
As $P_0 (x) = 1$, a value for $A_0$ can be found. It is
$$A_0 = \int^1_{-1} [P_0 (x)]^2 \, dx = \int^1_{-1} dx = 2.$$
So we have
\begin{align*}
A_n &= \frac{2n - 1}{2n + 1} A_{n - 1}\\
&= \frac{2n - 1}{2n + 1} \cdot \frac{2n - 3}{2n - 1} A_{n -
2}\\
& \hspace{1.0cm} \vdots\\
&= \frac{2n - 1}{2n + 1} \cdot \frac{2n - 3}{2n -
1} \cdots \frac{3}{5} \cdot \frac{1}{3} A_0\\
&= \frac{2n - 1}{2n + 1} \cdot \frac{2n - 3}{2n -
1} \cdots \frac{3}{5} \cdot \frac{1}{3} \cdot 2\\
&= \frac{2}{2n + 1}
\end{align*}
Hence
$$\int^1_{-1} [P_n (x)]^2 \, dx = \frac{2}{2n + 1},$$
and completes the proof.
Best Answer
I think the best way to approach this is as follows, note that \begin{align} (x^2 -1 ) = \frac{P_2 - 2}{3} \end{align} You can then use the following definition \begin{align} P_kP_l = \sum_{m=|k-l|}^{k+l} \begin{pmatrix} k & l & m \\ 0 & 0 & 0 \end{pmatrix}^2 (2m+1)P_m \end{align} This allows the integral to be written as follows \begin{align} \int_{-1}^{1} (x^2-1)^3P_iP_jP_k \; dx &= \int_{-1}^{1} \frac{1}{9}\left(P_2^3 + . . .-8 \right) P_i P_j P_k \; dx \end{align} The most difficult term to deal with is the $ P_2^3 P_i P_j P_k$ \begin{align} P_2^3 P_i P_j P_k &= \sum_{m=0}^{4} \begin{pmatrix} 2 & 2 & m \\ 0 & 0 & 0 \end{pmatrix}^2 (2m+1)P_m P_2 P_i P_j P_k \\ &= \sum_{m=0}^{4} \begin{pmatrix} 2 & 2 & m \\ 0 & 0 & 0 \end{pmatrix}^2 (2m+1) \sum_{n=|m-2|}^{m+2} \begin{pmatrix} 2 & m & n \\ 0 & 0 & 0 \end{pmatrix}^2 (2n+1)P_n P_i P_j P_k \\ &= \sum_{m=0}^{4} \begin{pmatrix} 2 & 2 & m \\ 0 & 0 & 0 \end{pmatrix}^2 (2m+1) \sum_{n=|m-2|}^{m+2} \begin{pmatrix} 2 & m & n \\ 0 & 0 & 0 \end{pmatrix}^2 (2n+1) \sum_{l=|n-i|}^{n+i} \begin{pmatrix} n & i & l \\ 0 & 0 & 0 \end{pmatrix}^2 (2l+1) P_l P_j P_k \end{align} Which can then make use of the usual triple integral. All other terms can be solved for in a similar manner.