Abstract: This study addresses the computational challenges associated with solving neutron transport equations in three-dimensional (3D) reactor cores, particularly in complex geometries where axial dimensions dominate radial dimensions. The method of characteristics (MOC) is widely used for their high precision, strong geometric adaptability, and suitability for parallel computing. However, conventional MOC faces significant limitations, including the requirement for fine axial grid discretization, which leads to high computational and storage demands. Additionally, while two-dimensional MOC coupled with axial methods reduces resource requirements, it introduces geometric restrictions and potential instabilities due to coupling. To overcome these limitations, this study aims to develop an advanced deterministic method that combines the advantages of MOC with reduced computational overhead. By focusing on high-order axial polynomial expansion and simplified radial modeling, the study seeks to enhance computational efficiency while maintaining accuracy, particularly for reactors with elongated fuel rods or plates. A novel characteristic plane method based on high-order axial polynomial expansion was developed. The method simplified radial expansion to zero order while enabling fourth-order axial polynomial expansion. This method leveraged the observation that in typical reactors, the axial dimensions of fuel rods or plates were much larger than their radial dimensions, making axial grid refinement more critical. The mathematical derivation involved two distinct scenarios of neutron transport within characteristic surfaces, deriving response matrices, and establishing flux representations. An iterative calculation framework was designed to handle boundary conditions and enable parallel computation. The method was implemented in a computational program to validate its effectiveness. Numerical validation using the Takeda and C5G7 benchmarks demonstrates the method’s accuracy and efficiency. In the Takeda benchmark, the axial grid size is successfully increased to 10 cm, matching the core height without additional layering, while maintaining acceptable deviations in k_eff (e.g., −38.3 pcm for case1 and −2.5 pcm for case2). Similarly, in the C5G7 benchmark, the axial grid size is extended to 14.28 cm, with k_eff deviations within acceptable ranges (e.g., −66.5 pcm for C5G7-3D and −112 pcm for Unrodded). The method significantly reduces computational time and storage requirements. The proposed high-order axial polynomial expansion for characteristic plane method offers a robust solution for 3D neutron transport calculations. By simplifying radial expansion and enhancing axial resolution, the method achieves grid sizes comparable to traditional block methods (e.g., 10-21.42 cm), significantly improving computational efficiency without compromising accuracy. This method is particularly advantageous for reactors with elongated axial dimensions, providing a practical alternative to conventional methods that require finer grids. The method’s stability, efficiency, and accuracy make it a promising tool for modern reactor design and analysis.