In a previous paper [Appl. Opt.52, 4151 (2013)], we presented the first- and second-order derivatives of a ray for a flat boundary surface to design prisms. In this paper, that scheme is extended to determine the Jacobian and Hessian matrices of a skew ray as it is reflected/refracted at a spherical boundary surface. The validity of the proposed approach as an analysis and design tool is demonstrated using an axis-symmetrical system for illustration purpose. It is found that these two matrices can provide the search direction used by existing gradient-based schemes to minimize the merit function during the optimization stage of the optical system design process. It is also possible to make the optical system designs more automatic, if the image defects can be extracted from the Jacobian and Hessian matrices of a skew ray.