A new numerical approach is attempted to solve the three-dimensional multigroup diffusion equation in rectangular geometry using the principles of variational calculus. The method essentially consists in separating the variables along the x, y, and z directions so that the neutron flux can be synthesized into three components to get an approximate solution. This solution can be used as a starting trial function for flux in the next step, where the variables are treated nonseparable in order to get a better solution. The effect of separation of variables is also studied. Detailed calculations are performed for a benchmark problem and the results are presented along with published values that are calculated using finite difference codes.