We solve the neutron diffusion equation by a wavelet Galerkin scheme in this paper. Wavelet functions are generated by dilation and translation operation on a scaling function. The wavelet functions are localized in space and have a recursive property, so these properties may be utilized to solve a differential equation that has severe “stiffness. ”The wavelet Galerkin method (WGM) represents the solution as a summation of Daubechies’ scaling functions, which are also used as the weighting function. The Daubechies’ scaling functions have the properties of orthogonality and high smoothness. Unlike the finite element method, the weighting function is the Daubechies’ scaling function, and the unknowns determined are not the fluxes of the nodes but the coefficients of the scaling functions. The scaling functions are overlapping in the nodes and require special treatment at interfaces between nodes and at the boundaries. We tested the WGM with several diffusion theory problems in reactor physics. The solutions are very accurate with increasing Daubechies’ order and dilation order. The boundary conditions are also satisfied very well. In particular, the WGM provides very accurate solutions for heterogeneous problems in which the flux distribution exhibits very steep gradients.We conclude that it is worthwhile investigating further the WGM for reactor physics problems and that numerical integration and acceleration of the matrix equation must be improved so as to reduce computing time.