In this paper a method is described that takes the nonlinear dynamic stiffness and damping coefficients for multiple hydrostatic bearings and incorporates them into a rotordynamic FEM model for a rotating machine. A Newton-Raphson iteration scheme is presented that uses updated bearing coefficients at every iteration to the solution. A non-linear computer program was written using the method described which models transient and synchronous response and calculates damped eigenvalues.