A perturbation method is developed for estimating the change in the solution of a reactive system to any order for a perturbation in the boundary condition in diffusion theory. The method derived gives formalisms for the eigenvalue, normalized flux, and homogenized parameters. Five examples are provided to verify the method as well as analyze the errors associated with it. The first example is very simple and solves the state of the system up to eighth order and gives a simple numerical analysis of a large perturbation. The next example gives an analytical solution up to second order. A two-region example is also given, which is partially numerical and partially analytical. An albedo test example shows that the higher-order terms all appear to be present in the formalism. The final example presents a simplified one-dimensional boiling water reactor core analyzed up to third order numerically. Applications of this method, error propagation, and future work are also discussed.