Кажется, ты застрял в классической петле, которая всё тормозит. Начни с того, чтобы векторизовал обновление всей сетки с помощью NumPy или CuPy, если есть GPU — это убьет накладные расходы Python. Потом добавь Numba JIT для самых внутренних циклов — пара декораторов, и увидишь ускорение в 5–10 раз. Если всё равно не поможет, попробуй переписать критически важный решатель на C++ и обернуть его с помощью pybind11. И не забудь заранее выделить все массивы и избегать динамического выделения памяти внутри ядра. Это должно дать толчок производительности.