Compared to traditional transverse integration methods, the variational nodal method, with its unique advantages, is more suitable for high-fidelity calculations of reactor physics in reactors with complex geometries and finer detail descriptions. In this study, the basic theory of the variational nodal method was derived and the VINUS code is developed. The neutron solver based on this method is adaptable to various geometric models, and showcased the code's fundamental framework. On this basis, a set of self-designed macroscopic cross-section benchmarks, actual macroscopic cross-section benchmark VVER-440, and few-group microscopic cross-section benchmark RBEC-M for fast spectrum reactors were used to verify different functionalities of VINUS. The results were shown that VINUS code maintains good computational accuracy and convergence trends. For the VVER-440 benchmarks, the deviation of keff of VINUS from reference is less than 100 pcm, and the maximum power deviation is less than 4 %. For the RBEC-M, the deviation of keff is 125 pcm, and the maximum power deviation is less than 5 %. These outcomes collectively demonstrate the solver's potential for engineering applications in future advanced reactor designs.