Неявная разностная схема ур-я движения в цилиндр. коорд.

[thunar]

Для декатовых координат неявная разностная схема решения уравнения движения заряженной частицы
$ \begin{array}{l} \frac{d^{2}\vec{r}}{dt^{2}}=\vec{E}+\vec{v}\times\vec{H}\\ \frac{d\vec{r}}{dt}=\vec{v} \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} $,
где индексом 0 отмечены значения с текущего временного слоя, а 1 -- со следующего. Имеем 6 линейных уравнений которые легко разрешается относительно $r_1$ и $v_1$. Получившаяся система решается с помощью нескольких последовательных итераций. На первой итерации значения полей с индексами 1 приравниваются к значениям полей с индексом 0, а затем берутся с прошлой итерации. Всё хорошо проблем нет.

Теперь, для цилиндрической системы координат с азимутальной симметрии уравнения движения (по координатам) записываются
$ \begin{array}{l} \frac{d^{2}R}{dt^{2}}=E_{r}+v_{\Theta}H_{z}+\frac{v_{\Theta}^{2}}{R}\\ \frac{d^{2}z}{dt^{2}}=E_{z}+v_{\Theta}H_{R}\\ \frac{d v_\Theta}{dt}=v_{z}H_{R}-v_{R}H_{z}-\frac{v_{\Theta}v_{R}}{R}\\ \frac{dR}{dt}=v_R\\ \frac{dz}{dt}=v_z \end{array} $,
где $R$ -- радиус, $z$ -- осевая координата, $v_\Theta$ -- азимутальная скорость.
Как в этом случае адекватно записать разностную схему и с учётом этих квадратичных членов разрешить её?