Calculating expectation value of the Hamiltonian squared 
So the main idea of the problem was to find the error in the argument, which I think I have a good grasp of. Basically, the Hamiltonian of the wavefunction is a constant non-zero value inside the box and zero outside. So $H^2\not=0$ at the edges, and hence the first statement saying that $H^2=0$, is definitely wrong. My question, is there a way to actually calculate the expectation value of $H^2$ because whenever I try to do that, I seem to end up with the derivative of a delta function and I have no idea how to figure that out.
Courtesy: Rudolf Ortvay Problems in Physics.
 A: The answer is that the considered function, though it belongs to the (essential)selfadjointness domain of $H$, 
$$D(H):= \{f: [-a,a] \to \mathbb C \:|\: f\in C^2([-a,a])\:, f(\pm a)=0\}$$
it does not belong to any (essential)selfadjoitness domain of $H^2$. So it does not make sense to compute its expectation value through that formula. To check my assertion try, integrating by parts, to prove that
$$\langle \Phi, H^2 \Psi\rangle=\langle H^2\Phi,  \Psi\rangle\qquad \Psi,\Phi\in D(H)\quad (false)$$
 You will see that the operator is not even symmetric on that domain because you can find functions vanishing at $\pm a$ though giving rise to non-vanishing boundary terms.
The same argument proves that $$0=\langle \Psi, H^2 \Psi\rangle\neq \langle H\Psi,  H\Psi\rangle>0$$
where $\Psi$ is that in OP's question.
I stress that there are however elements of $D(H)$ (a dense set actually!) which belong to an essentially selfadjointness domain of $D(H^2)$, these are in particular the finite linear combinations of  eigenvectors of $H$.   However $\Psi$ is not of that type.
This question concerns more mathematics than physics. It seems a bit, say, anomalous finding it in Rudolf Ortvay Problems in Physics.
A: Edit: The issue pointed out in the problem statement has implications in the below as well. Their failure might shed context on particular details of the problem. 
The second derivatives of the steady state solutions are zero at the boundary, yet the second derivative of the wave function is non-zero. What's going on here? 
Part 1: Assuming Hermiticity as in th below  has some problems.
$<\hat{H}^2>=\int_{-a}^a\Psi^*\hat{H}^2\Psi dx=\int_{-a}^a(\hat{H}\Psi)^*(\hat{H}\Psi)dx$
$\Psi(x)=N(a^2-x^2)$
$\hat{H}\Psi=\frac{-\hbar^2}{2m}(-2N)=\frac{h^2N}{m}$
$<\hat{H}^2>=\frac{\hbar^4N^2}{m^2}\int_{-a}^adx=\frac{2a\hbar^4N^2}{m^2}$
Specifically, things are thrown off by a boundary term:
$\Psi^*\frac{\partial^3 \Psi}{\partial x^3}-(\frac{\partial\Psi^*}{\partial x}\frac{\partial^2 \Psi}{\partial x^2})|_{-a}^a$
The wave function is an eigenvector of the second derivative so the boundary term can be rewritten as : $\sum_i[\Psi_i^*(\frac{-2mE_i}{\hbar^2})\frac{\partial \Psi_i}{\partial x}-(\frac{\partial \Psi^*}{\partial x})(\frac{-2mE_i}{\hbar^2})\Psi_i]=\frac{2m}{\hbar^2}\sum_i E_i(\Psi_i\frac{\partial \Psi^*_i}{\partial x}-\Psi_i^*\frac{\partial \Psi}{\partial x})$
Part 2: Plancherel's Theorem should break down as well.
$\Psi=\sum_{k=1}^\infty c_k \sqrt{\frac{1}{a}}\sin{(\frac{k\pi x}{2a})}$
$\hat{H}\Psi=\sum_{k=1}^\infty c_k\sqrt{\frac{1}{a}}\frac{\hbar^2k^2\pi^2}{8ma^2}\sin{(\frac{k\pi x}{L})}$
$<E>=\sum_{k=1}^\infty |c_k|^2| \frac{\hbar^2k^2\pi^2}{8ma^2}$
$<E^2>=\sum_{k=1}^\infty |c_k|^2| \frac{\hbar^4k^4\pi^4}{64m^2a^4}$
The orthonormality condition should yield the $c_k$'s:
$c_k=\int_{-a}^a N(a^2-x^2)\sqrt{\frac{1}{a}}\sin{(\frac{k \pi x}{2a})}dx$
Regardless, some conditions underlying the Plancherel Theorem are lacking casting doubts on the math. 
An important question to ask is, under what conditions does a series fail to represent a function? I suspect the  boundaries throw off the representation.
