  The non-linear transformations incurred by the rays in an optical system can
be suitably described by matrices to any desired order of approximation. In
systems composed of uniform refractive index elements, each individual ray
refraction or translation has an associated matrix and a succession of
transformations correspond to the product of the respective matrices. This
paper describes a general method to find the matrix coefficients for
translation and surface refraction irrespective of the surface shape or the
order of approximation. The choice of coordinates is unusual as the orientation
of the ray is characterised by the direction cosines, rather than slopes; this
is shown to greatly simplify and generalise coefficient calculation. Two
examples are shown in order to demonstrate the power of the method: The first
is the determination of seventh order coefficients for spherical surfaces and
the second is the determination of third order coefficients for a toroidal
surface.
