В статье предлагается численный алгоритм построения сплайн-функций Ляпунова для исследования абсолютной устойчивости нелинейных нестационарных систем. В случае асимптотической устойчивости системы выполнение алгоритма приведет к построению поверхности уровня функции Ляпунова в виде кусочно-гладкой (гладкой, при выполнении дополнительных условий) замкнутой поверхности размерности, равной размерности исходной системы. Показано, что данный алгоритм позволяет существенно улучшить оценки границы устойчивости, получаемые с помощью частотных методов. В отличие от алгоритма построения кусочно-линейных функций Ляпунова, время работы предлагаемого алгоритма построения сплайн-функций Ляпунова не стремится к бесконечности при приближении системы к границе устойчивости. Данное обстоятельство позволяет использовать модифицированный алгоритм для определения устойчивости систем, находящихся близко к границе устойчивости. Приведены оценки точности определения области устойчивости на примере конкретной системы 3-го порядка. Даны рекомендации по выбору начальных условий работы алгоритма. Ключевые слова: дифференциальные включения, нелинейные нестационарные системы, абсолютная устойчивость, функции Ляпунова, области устойчивости, сплайн Безье, полиномы Бернштейна.