In this work the paraxial optical imaging is generally described by means of three square matrices: one unitary system matrix and two operation matrices with determinants equal the Lagrange-Helmholtz invariant. Elements of system matrix are functions of design parameters while elements of the operation matrix depend on input and output coordinates of characteristic rays. Each matrix has only three independent elements. Internal system parameters are determined from equations created of system matrix elements, which values are dependent on the operation matrix. Matrix approach enables the solution of only three non-linear equations with respect to system parameters. Matrix approach has also another advantage. It enables the determination of number of degrees of freedom. We have a superiority of parameters over the number of equations when the number of components is bigger than 2. The more complex is model the higher degree of freedom it has. There are special ways of reducing the number of degree of freedom: by selection of spaces between component, introduction of additional requirements and criteria of distribution of optical powers. Significant help is in defining all the spaces between components, what means full control of the components position and their dimensions. In such a case the only thing left is the determination of optical powers, while the number of degree of freedom is equal k-1 (k is the number of components). In this work the computer program realizing described algorithms has been developed. This program was tested with specially selected examples. Results of calculation for two interesting applications are also given.