A non-linear finite element calculation method based on the coupling and nonlinearity of bearing stiffness is presented to solve the difficulties for the analysis of multi-support shaft system. A bearing element is built in the present method. The coupling and nonlinearity of bearing stiffness, the preload of bearing, the roller taper angle, and the other details are considered, and the analytical expression of bearing stiffness matrix is deduced. The shaft is modeled by Euler beam element with shear deformation, and system computational model is assembled by the model of shaft and bearings. The NR method is used to solve the model. An example shows that this method can calculate the load of each bearing in every direction and the shaft deformation for multi-support shaft system by given external load correctly. This method has high computational efficiency, and is suitable for the analysis of multi-support shaft system.