Linear vibration analysis has been widely used in various industries such as aero, auto, and electronic industries for decades, and we have two challenges. One is the increasing demand for higher frequency ranges to solve. For example, Noise, Vibration, and Harshness (NVH) for eMobility require considering electric motor vibration, road noise, and wind noise up to 1kHz or more. Another is to take into account nonlinearity in the vibration analysis. Practically, any applicative structures include nonlinearity, such as bushings, bearings, squeeze film dampers, and contact. Traditionally, the real mode-based reduction is applied to the equation of motion to accelerate the solution of the equation. However, it will produce a full matrix if the damping, gyro effect, or fluid-structure coupling is considered which leads to slow processing in the response analysis with many real modes. To address this issue, a new type of complex mode-based reduction method has been developed by Kyushu University, which produces a completely diagonal matrix of the equation of motion. The solution of the reduced equation is drastically faster than the traditional real mode-based reduced equation. Higher order modes are required to solve nonlinear harmonic responses than to solve linear harmonic responses. To reduce the number of complex modes to be extracted, the correction terms of the truncated modes are considered. Moreover, the part of complex modes that do not affect the response much can be deleted from the equation to speed up the solution more with good accuracy by adding the correction term, which is calculated using the complex modes, to the nonlinear degree of freedom. The new type of complex mode-based reduction method and the correction terms are introduced in a toolkit working with MSC Nastran to solve for the nonlinear harmonic response. The nonlinear harmonic response is solved by the shooting method which is the method for a boundary value problem by reducing it to an initial value problem. The toolkit is applied to a structure with nonlinear bushings and a structure with contacting another structure. The harmonic responses are solved successfully, and the characteristics of the harmonic response are classified as stable or unstable. The robustness and performance of the method will be demonstrated with aerospace and automotive examples.