Уравнение лапласа в разностном представлении

[binSmile]

Добрый день!
Проблема, не могу корректно преобразовать уравнение Лапласа в цилиндрической системе координат в разностное представление, для решение методом сеток. Задача у меня осесимметричная и следовательно я работаю с уравнением Лапласа в двух координатах: r и z.
Изначально, уравнение выглядит так:
$\frac{1}{r}\frac{\partial }{\partial r}(\varepsilon r \frac{\partial U}{\partial r}) + \frac{\partial }{\partial z}(\varepsilon \frac{\partial U}{\partial z})=0$

Сетку я строю так, что изменение узла по i - это изменение по оси r, а j - по оси z.
$\varepsilon$ - у меня не константа.
Я следовал по такой логике, что первая производная, это разность соседних точек. следовательно производная от производной - разность двух разностей %)
Соответственно, если не выносить общие множители, у меня получается так:


$\frac{1}{r}(\frac{1}{h_r}(\frac{\varepsilon_{i+1}r_{i+1}(U_{i+1}-U)}{h_r}-\frac{\varepsilon r(U-U_{i-1})}{h_r}))+\frac{1}{h_z}(\frac{\varepsilon_{j+1}(U_{j+1}-U)}{h_z}-\frac{\varepsilon (U-U_{j-1})}{h_z}) = 0$

где $U=U_{i,j}, U_{i+1} = U_{i+1,j}$ аналогично для индекса $j$ и функции $\varepsilon $

Тоже, если вынести общие множители U:
$U_{i+1,j} \frac {\varepsilon_{i+1,j}(i+1)}{h^2_r i} + U_{i-1}\frac{\varepsilon_{i,j}}{h^2_r}+U_{i,j+1} \frac{\varepsilon_{i,j+1}}{h^2_z} + U_{i,j-1}\frac{{\varepsilon_{i,j}}}{h^2_z}-U(\frac{\varepsilon_{i+1,j}(i+1)}{h^2_ri}+\frac{\varepsilon_{i,j}}{h^2_r}+\frac{\varepsilon_{i,j+1}}{h^2_z}+\frac{\varepsilon_{i,j}}{h^2_z})=0$

Собственно, возможно, я ошибаюсь в дальнейших вычислениях. Но проверить не закралась ли тут ошибка, я не придумал варианта как.

Заранее благодарен за советы.