Two-point boundary-value ray tracing in anisotropic elastic media, based on the ray bending method, is a highly non-linear problem. It can be solved by the Newton method, which requires first and second spatial, directional and mixed derivatives of the ray (group) velocity at each node along a trial ray trajectory between two fixed endpoints. The second derivatives also provide the curvature components of the propagating wave fronts and thus make it possible to compute dynamic characteristics (e.g. geometrical spreading) along the rays, as well as paraxial rays and beams. We have developed an original novel methodology for obtaining analytically these spatial/directional gradient vectors and spatial/directional/mixed Hessian matrices for general anisotropic (triclinic) media. The derivatives for higher crystal' symmetry classes, such as transverse isotropy, orthorhombic and monoclinic media (including those with tilted symmetry axes or planes) are explicitly provided; each is followed by a numerical example. We validate our exact analytic derivatives by comparing them to numerical approximations computed with finite difference schemes.