Power system (PS) stability analysis is one of the fundamental criteria in system operation analysis. The issue of power system (PS) stability is one of the fundamental issues that is addressed by appropriate system stability analysis. Power system stability simulation programs in the time domain are an integral part of the power system analysis tools. With time domain simulation programs, the dynamic response of the power system to disturbances or changes in the system state can be calculated using appropriate mathematical models based on numerical algorithms. The developed numerical model for the analysis of angular stability in the power system (PS) is based on the finite element technique and time-varying phasors. The complex numerical analysis was simultaneously carried out in the time and phasor domains. The basis of the developed numerical model for the analysis of angular stability in the PS is the application of the finite element technique to the PS so that the considered system is divided into smaller parts, which are treated in the calculation process as separate finite elements with a certain number of local nodes. As the most important part of the power system, in the first step, a subtransient numerical model of a synchronous generator was developed, which is defined and treated as a finite element with three local nodes. The developed numerical model of a synchronous generator for the analysis of angular stability in the power system also includes the originally developed numerical model of the turbine regulator and the originally developed numerical model of the excitation regulator. Also, based on the appropriate mathematical models, numerical models of other parts of the power system (transformers, lines, equivalent network) were developed, which were treated as separate finite elements in the calculation process. The system of differential equations of the generator, turbine regulator and excitation regulator is included in the numerical model by using the numerical integration of the aforementioned differential equations using the Heun method and the generalized trapezoidal rule (ϑ - method). Therefore, with the help of this procedure, i.e. the Heun and ϑ - methods, the system of differential equations is reduced to a system of algebraic equations in one time increment. The accuracy of the developed numerical model for the analysis of angular stability in the power system was confirmed by comparing the calculation results with the results obtained using the EMTP (Electromagnetic Transient Program) software package. The calculation results obtained using the Heun method are more accurate than the results obtained using the ϑ - method.