The two most common optical boundaries in geometrical optics are the spherical and flat. The present group recently derived the third-order derivative matrix of a skew ray with respect to the source ray vector for a ray reflected/refracted at a flat boundary. The proposed method was based on a differential geometry approach, and hence had the advantages of an improved accuracy and the need to trace just one ray. In the present study, the method is extended to a ray incident on a spherical boundary. The derived matrix is used to explore the primary wavefront spherical aberration of an axis-symmetrical system. Its result is identical to that obtained from Zemax simulations. The estimation capability of the fifth-order Taylor series expansion of a ray is investigated by using the finite difference methods and the developed matrix. The proposed matrix can serve as a useful basis for determining the higher-order differential derivative matrices of a ray to explore higher-order ray or wavefront aberrations in future studies.