A method for determining the mathematical adjoint solution of a higher order nodal expansion method (NEM) based on the simultaneous solution of multigroup equations for each node in the rectangular geometry is presented. In the higher order NEM, the forward NEM equations in a given node include not only the nodal balance and interface-current equations but also weighted residual method (WRM) equations for higher order expansion coefficients. In deriving the mathematical adjoint equations corresponding to these forward NEM equations, the transverse leakage terms in the WRM equations need to be replaced by partial currents. Because transverse leakage terms of a node are linked to partial currents of many neighboring nodes, replacement of transverse leakage terms by partial currents results in complicated WRM equations. Because mathematical adjoint equations are obtained by transposing the nodal forward equations, direct use of these complicated WRM equations makes the numerical computation of the adjoint solution inefficient. This problem is avoided by treating the transverse leakage terms contained in the WRM equations as additional unknowns and by including the equations defining the transverse leakage terms in terms of partial currents into the nodal forward equations. The mathematical adjoint equations are then derived by transposing the resulting nodal forward equations. This adjoint solution method is verified by comparing nodal adjoint fluxes with the fine-mesh VENTURE solution for the International Atomic Energy Agency (IAEA) pressurized water reactor (PWR) benchmark problem and by comparing the local reactivity changes computed with first-order perturbation theory for the IAEA PWR and the Yonggwang unit 2 PWR with the exact reactivity values determined from the eigenvalue difference between perturbed and unperturbed cores.