This work develops efficient algorithms for numerically solving the neutron diffusion equation by a wavelet Galerkin method (WGM). One of the main problems encountered in solving neutron diffusion equation using WGM is the treatment of the boundary and interface conditions. In one-dimensional problems, the boundaries of the wavelet series expansions are assumed to be the analytical boundaries of the problem, and the boundary condition equations are replaced by end equations in Galerkin system. In two-dimensional problems, in order to maintain the integrity of the system, the boundaries of the wavelet series are shifted until the end is independent on any expansion coefficients of the scaling function that affect the solution within the real boundaries.Since the scaling functions are compactly supported, only a finite number of the connection coefficients are nonzero. The resultant matrix has a block diagonal structure, which can be inverted easily, therefore, enables us to extend the solution to two-dimensional heterogeneous cases and make the inner iteration efficient in the eigenvalue and multigroup problems. We tested our WGM algorithm with several diffusion problems including two-dimensional heterogeneous problems and multi-group problems. The solutions are very accurate with a proper selection of Daubechies' order and dilation order. Our WGM algorithm provides very accurate solutions for one-dimensional and two-dimensional heterogeneous problems in which the flux exhibits very steep gradients.