A method is described that allows the computation and analysis of high-order spin maps for general non-autonomous optical systems. It is shown how the equations of motion in curvilinear coordinates resulting from the Thomas-BMT equation can be solved within a differential algebraic framework in SU(2) and SO(3) representations. The resulting maps are subjected to a spin-orbit normal form transformation, and the nonlinear orbit dependencies of the invariant polarization axis as well as the orbit dependent spin tune can be obtained. For the case of electron machines, the resulting invariant polarization can be used to determine the radiative equilibrium polarization via the Derbenev-Kondratenko approach. Both the computation of the spin map as well as the algorithm for the computation of the invariant axis have been implemented in the code COSY INFINITY [2] [3] [4].