A first-principles method is presented to calculate elastic constants up to the fourth order of crystals with the cubic and hexagonal symmetries. The method relies on the numerical differentiation of the second Piola-Kirchhoff stress tensor and a density functional theory approach to compute the Cauchy stress tensors for a minimal list of strained configurations of a reference state. The number of strained configurations required to calculate the independent elastic constants of the second, third, and fourth order is 24 and 37 for crystals with the cubic and hexagonal symmetries, respectively. Here, this method is applied to five crystalline materials with the cubic symmetry (diamond, silicon, aluminum, silver, and gold) and two metals with the hexagonal close packing structure (beryllium and magnesium). Our results are compared to available experimental data and previous computational studies. Calculated linear and nonlinear elastic constants are also used, within a nonlinear elasticity treatment of a material, to predict values of volume and bulk modulus at zero temperature over an interval of pressures. To further validate our method, these predictions are compared to results obtained from explicit density functional theory calculations.