  The model and algorithm for the cooling of the magnetized neutron stars are
presented. The cooling evolution described by system of parabolic partial
differential equations with non-linear coefficients is solved using Alternating
Direction Implicit method. The difference scheme and the preliminary results of
simulations are presented.
