Для декатовых координат неявная разностная схема решения уравнения движения заряженной частицы

записывается как
![$
\begin{array}{l}
\vec{v}_{1}=\vec{v}_{0}+0.5\delta t\left(\vec{E}_{1}+\vec{E}_{0}+\left[\vec{v}_{0}\times\vec{H}_{0}\right]+\left[\vec{v}_{1}\times\vec{H}_{1}\right]\right)\\
\vec{r}_{1}=r_{0}+0.5\delta t\left(\vec{v}_{0}+\vec{v}_{1}\right)
\end{array}
$ $
\begin{array}{l}
\vec{v}_{1}=\vec{v}_{0}+0.5\delta t\left(\vec{E}_{1}+\vec{E}_{0}+\left[\vec{v}_{0}\times\vec{H}_{0}\right]+\left[\vec{v}_{1}\times\vec{H}_{1}\right]\right)\\
\vec{r}_{1}=r_{0}+0.5\delta t\left(\vec{v}_{0}+\vec{v}_{1}\right)
\end{array}
$](https://dxdy-03.korotkov.co.uk/f/a/0/3/a03a780444afc29a8890609e23b5373682.png)
,
где индексом 0 отмечены значения с текущего временного слоя, а 1 -- со следующего. Имеем 6 линейных уравнений которые легко разрешается относительно

и

. Получившаяся система решается с помощью нескольких последовательных итераций. На первой итерации значения полей с индексами 1 приравниваются к значениям полей с индексом 0, а затем берутся с прошлой итерации. Всё хорошо проблем нет.
Теперь, для цилиндрической системы координат с азимутальной симметрии уравнения движения (по координатам) записываются

,
где

-- радиус,

-- осевая координата,

-- азимутальная скорость.
Как в этом случае адекватно записать разностную схему и с учётом этих квадратичных членов разрешить её?