Dynamics of multi-cylinder engines have been performed by the conventional methods using lumped pure torsional systems. This article offers a finite element method of performing elastodynamic analysis of multi-cylinder engines considering them as spatial linkage systems with true spatial geometries of the crankshaft and the linkage loops. An engine can have any number of cylinders with linear offsets and angular orientations relative to each other. A three-dimensional finite-line element with isoparametric joint irregularity freedoms is developed and used. Consistent or lumped mass systems can be used. Elastodynamics of engines is considered in two forms: (a) kinetoelastostatics (KES) where all forces and torques acting on the system are considered except the vibratory motion effect; (b) kinetoelastodynamics (KED) where the forced and damped equations of motion of the system are solved. Matrix exponential method of solving equations of KED motions are presented and used. It is proven to be a very efficient and stable technique for the solutions of large systems of linear and nonlinear differential equations of any order. After solving for the generalized coordinates, time histories of the neutral coordinate displacements, forces, moments, stresses, bearing forces, and generated torque are determined for as many work cycles as desired. A generalized computer program performing KED and KES studies of any multi-cylinder engine is made available for industrial use. KED and KES analyses of a four-cylinder automobile engine are performed.