A methodology is presented that allows a higher-order accurate treatment of system perturbations that are assumed to have a substantial magnitude and therefore also a substantial effect on flux distributions in nuclear systems. Examples are localized material choice variations, burnable poison density variations at lattice level, complete fuel assembly permutations at core level, or specific uncertainties defined in the system composition. For these cases, it is necessary to raise the accuracy of the estimated responses above what can be achieved using first-order perturbation methods only, of course preferably without having to simply pursue computationally expensive exact recalculations for each case if the effects of many variations or uncertainties are to be assessed. Focusing on the neutronics of multiplying systems (without thermal-hydraulic feedback mechanisms incorporated), the setup of a polynomial form for quantification of the flux shape change due to imposed system perturbations is pursued. In a mathematical sense, this method allows one to set up a polynomial expansion of the change in the lowest-mode solution of the neutronics eigensystem due to an imposed perturbation in the operators determining the lowest-mode solution and eigenvalue. This form features the property that the flux shape change, caused by variations in certain parameters localized in space and energy, can be expanded polynomially up to higher-order accuracy, with the imposed system composition variations themselves as functional arguments. Numerical results, showing the validity of the method, are reported, and possible application areas are discussed.