A key difficulty in nuclear reactor analysis is to provide a wide description of nuclear systems through the coupling of neutronics and thermal hydraulics. COMSOL Multiphysics software provides an advanced tool to integrate user-defined physics modules to analyze reactors at steady and transient state by mathematics module. To validate the capability of solving different nuclear reactor phenomena e.g., neutron flux distribution and effective multiplication factor (keff ) in COMSOL software, it’s worth solving a benchmark problem. A classic 3D IAEA benchmark problem is selected for the estimation of effective multiplication factor, neutron flux, radial, and axial power distributions in the core using the neutron diffusion equation. The 3D model was developed with SolidWorks and afterward imported into COMSOL to solve two-group neutron diffusion equations. To avoid the peculiarities that commonly occur when solving the neutron diffusion equation with unstructured grids, an adaptive meshing approach was utilized. The effective multiplication factor, thermal and fast neutron flux profiles, as well as power distributions, were calculated and compared with the results of another standard PARCS code. There is good agreement between the COMSOL and PARCS code results. The value of keff is 1.02799 which shows a smaller difference of 0.11% compared to the reference value from the PARCS code. The results presented in this paper are primarily intended as a demonstration of the neutronics behavior of this benchmark problem with COMSOL multiphysics software. Based on the present study, it is evident that the COMSOL software could be used for neutronics analysis for nuclear reactors like PWR. VVER etc.