A General Method for the Determination of Matrix Coefficients for High   Order Optical System Modelling
  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. 
Jose B. Almeida