I propose a method for derivative approximation, virtually unknown in the geophysical literature but is rapidly gaining recognition in other areas of computational science. The technique, called the complex-step derivative (CSD), is based on the theory of complex variables and approximates first-order derivatives of real-valued functions that are analytic in the complex plane. I extend the CSD technique to second-order derivatives, while preserving the robustness of the first-order formula. Thus the formulas together provide a complete differentiation system (referred to as semiautomatic differentiation, or SD) that allows efficient and accurate approximation of gradients, Jacobians and Hessians, as well as 2D and 3D partial and cross-partial spatial derivatives. Performance evaluation tests indicate that in comparison with ordinary finite-difference schemes (FD), the SD scheme is six to eight orders of magnitude more accurate, numerically highly stable, and step-size insensitive, which are major advantages over FD. The method shares with FD its attractive feature of ease of implementation. The SD scheme is implemented in a MATLAB toolkit. The validity of the CSD method depends critically on the requirement that the target function of differentiation be analytic in the complex plane. Therefore, prior to using the CSD method, one must ensure that all functions in the complex library are defined so that they satisfy this requirement. Specialized complex-function libraries that resolve this and other technical issues for CSD applications are available in the public domain.