The use of a numerical differentiation formula (NDF) is an excellent method for solving stiff ordinary differential equations. However, the NDF method cannot fully adapt to all stiff systems. An optimal general method for optimizing NDF coefficients using a back-propagation neural network is proposed in this work that can be used for different systems of stiff equations. The ranges of stability of the first- to fourth-order coefficients are obtained by analyzing the definition of stability. In order to solve the different stiff systems by changing the NDF coefficients, the relationship between the eigenvalues of the Jacobian matrix and NDF coefficients is analyzed. The back-propagation neural network is used to describe the relationship between them and predict the optimal parameters of different stiff systems. Compared with the NDF coefficients, the absolute and mean square errors of the numerical solutions are smaller. The simulation results show that the numerical solution’s accuracy is higher. Because optimizing the NDF coefficients improves the accuracy of the simulation results, the new approach is more suitable for solving stiff ordinary differential equations.