стационарное ур-е непрерывности + ур-е Пуассона

[vJB]

Товарищи,

помогите идейно с задачей.
Физическая часть:
плазма находится в эл.ст. поле (геометрия - цилиндрическая), нужно определить:
- профили концентраций электронов $n_e(z,r)$, ионов $n_i(z,r)$ (азимутальная симметрия)
- распределение потенциала $\varphi(z,r)$

Математическая часть:
В общем случае система уравнений такова

$$ \left \{ \begin{array}{ll} \frac {\partial n_e}{\partial t} + \vec{ \nabla } \cdot \vec{ \Gamma }_e = 0 & (1) \\ \frac {\partial n_i}{\partial t} + \vec{ \nabla } \cdot \vec{ \Gamma }_i = 0 & (2) \\ \Delta \varphi = - \frac {e} {\varepsilon _0} \cdot \left ( n_i - n_e \right ) & (3) \end{array} \right. $$

условия на границе:
$ n_{e} (0,r) = n_{e0} $
$ n_{i} (0,r) = n_{i0} $
значение потенциала \varphi на границах задано

Потоки частиц

% electron flux $$ \vec{\Gamma}_e = - \mu_e \cdot n_e \cdot \vec{E} - D_e \cdot \vec{\nabla} n_e \eqno (4) $$ % ion flux $$ \vec{\Gamma}_i = \mu_i \cdot n_i \cdot \vec{E} - D_i \cdot \vec{\nabla} n_i \eqno (5) $$

Напряженность поля

$$ \vec{\nabla} \cdot \vec{E} = \frac{e}{\epsilon _0} \cdot \left ( n_i - n_e \right ) \eqno (6) $$

$$ (3) \Leftrightarrow (6) $$

В стационарном случае

$$ \frac { \partial n_e }{ \partial t } = 0, \quad \frac { \partial n_i }{ \partial t } = 0 $$

Система ур-й, которую нужно решить численно

$$ \left \{ \begin{array}{ll} \mu_e \cdot \left ( \left ( \vec{\nabla} n_e \right ) \cdot \vec{E} + n_e \cdot \left ( \vec{\nabla} \cdot \vec{E} \right ) \right ) + D_e \cdot \Delta n_e = 0 & (7) \\ \mu_i \cdot \left ( \left ( \vec{\nabla} n_i \right ) \cdot \vec{E} + n_i \cdot \left ( \vec{\nabla} \cdot \vec{E} \right ) \right ) - D_i \cdot \Delta n_i = 0 & (8) \\ \Delta \varphi = - \frac {e} {\epsilon _0} \cdot \left ( n_i - n_e \right ) & (3) \end{array} \right. $$

Подставляем $(6)$ в $(7),(8)$
Получим окончательно

$$ \left \{ \begin{array}{ll} \left ( \vec{\nabla} n_e \right ) \cdot \left ( - \vec{\nabla} \varphi \right ) + \frac{e}{\epsilon_0} \cdot n_e \cdot \left ( n_i - n_e \right ) + \frac{D_e}{\mu_e} \cdot \Delta n_e = 0 & (9) \\ \left ( \vec{\nabla} n_i \right ) \cdot \left ( - \vec{\nabla} \varphi \right ) + \frac{e}{\epsilon_0} \cdot n_i \cdot \left ( n_i - n_e \right ) - \frac{D_i}{\mu_i} \cdot \Delta n_i = 0 & (10) \\ \Delta \varphi = - \frac {e} {\epsilon _0} \cdot \left ( n_i - n_e \right ) & (3) \end{array} \right. $$

Получили систему нелинейных ур-й в ЧП.
Далее, я записываю эти ур-я в конечных разностях (в цилиндрических координатах) и решаю методом Бройдена.
Но - решение расходится (норма вектора решения возрастает с каждой итерацией до переполнения).

Вопрос:
Кто-нибудь сталкивался с такими задачами? Как решали?
Может, хоть ссылками на соответствующие статьи/литературу поможете?
p.s.
картинка с хостинга не вставляется: "Не удалось определить размеры изображения"