In this work, we present an algorithm that enables computation of inertial modes and their corresponding frequencies in a rotating triaxial ellipsoid. The method consists of projecting the inertial mode equation onto finite-dimensional bases of polynomial vector fields. It is shown that this leads to a well-posed eigenvalue problem, and hence, that eigenmodes are of polynomial form. Furthermore, these results shed new light onto the question whether the eigenmodes form a complete basis, i.e. whether any arbitrary velocity field can be expanded in a sum of inertial modes. Finally, we prove that two intriguing integral properties of inertial modes in rotating spheres and spheroids also extend to triaxial ellipsoids.