First things first: the operator which corresponds to the energy is the Hamiltonian, typically written as $H$. So when you want to get the expectation value of the energy, you evaluate $\langle H\rangle$.
Now, there are multiple ways to do this. One way is to use the Schroedinger equation to get
$$\langle H\rangle = \left\langle i\hbar\frac{\partial}{\partial t}\right\rangle = \int\Psi^*(x,t) i\hbar\frac{\partial}{\partial t}\Psi(x,t)\,\mathrm{d}x\tag{1}$$
That calculation is completely general, i.e. it is valid in any situation for which the Schroedinger equation applies.
Another way to get $\langle H\rangle$ is to use the definition of the Hamiltonian operator, which in nonrelativistic QM is
$$H = \frac{p^2}{2m} + V(x,t) = -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x,t)$$
That gives you
$$\begin{align}\langle H\rangle &= \biggl\langle -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x,t)\biggr\rangle \\ &= \int\Psi^*(x,t)\biggl(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x,t)\biggr)\Psi(x,t)\,\mathrm{d}x\tag{2}\end{align}$$
Either (1) or (2) works in general.
For a free particle only, the potential $V(x,t)$ is zero, and you get
$$\langle H\rangle = \int\Psi^*(x,t)\biggl(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}\biggr)\Psi(x,t)\,\mathrm{d}x$$
which is the expression you saw on that web page.
In your problem, you need integrals of kind :
$I_{2n} = \int x^{2n} e^{- \large \frac{x^2}{a}} ~ dx$
Note first that $I_0 = (\pi)^\frac{1}{2} (\frac{1}{a})^ {-\frac{1}{2}}$
Now, it is easy to see that there is a reccurence relation between the integrals :
$$I_{2n+2} = - \frac{\partial I_{2n}}{\partial (\frac{1}{a}) } $$
For instance,
$$I_2 = - \frac{\partial I_{0}}{\partial (\frac{1}{a}) } = \frac{1}{2}(\pi)^\frac{1}{2} (\frac{1}{a})^ {- \large\frac{3}{2}} = \frac{1}{2}(\pi)^\frac{1}{2} ~a^ {\large\frac{3}{2}}$$
$$I_4 = - \frac{\partial I_{2}}{\partial (\frac{1}{a}) } = \frac{3}{2} \frac{1}{2}(\pi)^\frac{1}{2} (\frac{1}{a})^ {- \large\frac{5}{2}} = \frac{3}{2} \frac{1}{2}(\pi)^\frac{1}{2} ~a^ {\large\frac{5}{2}}$$
A general formula is :
$$I_{2n} = I_0 ~(2n-1)!! ~(\frac{a}{2})^n = \frac{(\pi)^\frac{1}{2}}{2^n} ~(2n-1)!! ~a^{n+\frac{1}{2}}$$
where $(2n-1)!! = (2n-1)(2n-3)......5.3.1$
Best Answer
Well, $\widehat{p^2} = \hat{p}^2= \hat{p} \hat{p}$.
So, in the position basis it is $-\hbar^2 \frac{d^2}{dx^2}$, and $\langle p^2 \rangle = \int_{-\infty}^\infty \bar{\Psi}\left(-\hbar^2 \frac{d^2}{dx^2} \right)\Psi dx$.
Note: $\hat{p}$ is technically not equal to $-i\hbar d/dx$, but rather in the position basis $\langle x | \hat{p}| x' \rangle = -i\hbar d/dx \delta(x-x')$.