Вообще-то универсальные алгоритмы для эволюционных систем уравнении существуют. Но у них есть проблема очень простая: чтобы они поработали, нужны уж очень маленкие степени времени.
Скажем,

. Там очень просто и универсально можно написать

и все это решается просто, все

вычисляются просто если все начальные данные

даны. Но есть проблема - маленькие колебания начальных данных могут вести к большим осциллациям. В принципе надо просто взять

достаточно маленьким. Но в практике слишком маленький

значит слишком много шагов, и поэтому нужно применить другие методы, которые менее универсальные, но при которых даже при

побольше не появятся никаких неустойчивостей.