Для моделирования тепловой эволюции нейтронной звезды построена численная схема, свободная от традиционного разделения звезды на сильно вырожденную внутреннюю часть и квазистационарную теплоизолирующую оболочку. Учтён ряд недавних уточнений микрофизики. Для нейтронных звёзд с сильными магнитными полями учтено влияние квантования Ландау на термодинамические функции, теплопроводность и структуру внешней коры. Выполнено численное моделирование остывания обычных нейтронных звёзд и магнитаров с неаккрецированными и аккрецированными оболочками. Проведено сравнение результатов моделирования и наблюдений. Изучена чувствительность остывания нейтронных звёзд к модели и толщине квазистационарной оболочки, к различным микрофизическим факторам, а также к эффектам, вызванным сверхсильными магнитными полями. Работа поддержана грантом РНФ 14-12-00316.