<p>Periodic responses of nonlinear dynamical systems are of prime importance for engineering vibration analysis. Their time-domain trajectories can be analytically represented via the synthesis of Fourier-based harmonics with undetermined coefficients. The paper aims to overcome two critical difficulties in solving periodic responses: the high computational cost of solving high-order harmonic coefficients for high-precision responses, and the difficulty in one-shot global capture of all coexisting stable and unstable periodic solutions. To this end, this paper presents a novel neural network-based subspace harmonic expansion (NNSHE) method for efficiently capturing all potential high-precision periodic responses. The method integrates Fourier harmonic synthesis into a physically constrained neural network framework, equivalently transforming the solution of undetermined harmonic coefficients into the iterative optimization of network hyperparameters. The weights space with Fourier harmonic bases is defined and discretized into cells. Utilizing the error tolerance of cells, a structurally extendable neural network framework is established to find harmonic coefficients from their low-dimensional and coarse covering set to higher dimensional space. Cell information provides effective prior guidance to accelerate weight optimization significantly. Meanwhile, a global penalty function is embedded into the loss function to avoid local minima and ensure full exploration of the cell space, enabling all coexisting periodic solutions to be captured in a single optimization. The performances of the proposed method in global response prediction and high-order (e.g., 500 orders) harmonic expansion are demonstrated through the analysis of two typical nonlinear systems.</p>