A number of improvements have been made to the Hermite method in order to obtain a high order finite element method capable of solving the neutron diffusion equation. First, a variational formulation of the equation is used to obtain a Weierstrass-Erdmann-type coupling relation valid at all points in the domain, singular and nonsingular. The basic solution yielded by this type of discretization is obtained by the inverse power method with variational acceleration of outer iterations. The linear systems appearing in the inverse power method are solved using a one-way dissection algorithm followed by asymmetric block factorization. These procedures were programmed in the BIVAC code for a treatment of the neutron diffusion equation with a two-dimensional reactor representation. The Hermite method was then compared with alternative approaches to a solution. The tests correspond to two-dimensional configurations of pressurized water reactors and CANDU reactors.