An automated procedure for calculating second-order elastic constants for crystalline systems of any symmetry using the CRYSTAL program is described. Second derivatives with respect to strain are evaluated numerically from analytical gradients. The internal co-ordinates are re-optimized with each applied strain. Point group symmetry is exploited to reduce the number of needed deformations according to Laue classes. A set of test cases covering many of the crystal classes is used to document the numerical accuracy of the scheme, and to define default values of the computational parameters so as to reduce the input file to a single keyword.