A new symbolic differentiation algorithm is proposed in this paper to automatically generate the inverse dynamics of flexible-joint robots in symbolic form, and results obtained can be used in real-time applications. The proposed method with computational complexity is developed based on the recursive Newton–Euler algorithm, the chain rule of differentiation, and the computer algebra system. The input of the proposed algorithm consists of symbolic matrices describing the kinematic and dynamic parameters of the robot. The output is the inverse dynamics solution written in portable and optimized code (C-code/Matlab-code). An exemplary, numerical simulation for inverse dynamics of the Kuka LWR4 robot with seven flexible joints is conducted using matlab, in which the computational time per cycle of inverse dynamics is about 0.02 ms. The numerical example provides very good matching results versus existing methods, while requiring much less computation time and complexity.