Как установить граничные условия для трехмерного уравнения Пуассона в Matlab

Question

Я пытаюсь вычислить электрический c потенциал в точке (x, y, z) путем решения приведенного ниже трехмерного уравнения Пуассона с использованием метода конечных разностей.

\frac{\delta^2 V}{\delta^2 x^2} +\frac{\delta^2 V}{\delta^2 y^2} +\frac{\delta^2 V}{\delta^2 z^2} 

У меня есть плотности заряда в различных положениях rho (x, y, z).

Ниже приведены граничные условия;

1. Граничное условие Дирихле применяется вверху и внизу плоскостей прямоугольной сетки angular.
2. Электричество c Потенциал должен быть включен установкой V (x , y, 0) = 0 и V (x, y, h) = -Eh, где h - высота блока моделирования.
3. Граничные условия Неймана также применяются на оставшихся интерфейсах блока путем установки в \frac{\delta V}{\delta x} грани с константой x, грани с константой \frac{\delta V}{\delta y} y и \frac{\delta V}{\delta z} грани с константой z.

Я сделал реализацию граничных условий с этим кодом, и я бы хотите убедиться в правильности реализации граничных условий.

y1 = linspace(0,10,20);
z1 = linspace(0,10,20);
V = zeros(length(x1),length(y1),length(z1)); 
%Dirichlet Boundary Conditions
%Top plane
V(:,end,end) = 0;
V(end,:,end) = 0;
V(end,end,:) = 0;
% Bottom plane
V(:,1,1) = 0;
V(1,:,1) = 0;
V(1,1,:) = 0;
%Incoporated electric potential
V(:,:,1)  = 0; 
V(:,:,end) = -40*z1(end);
%Neumann Boundary Condition
i  = 2:length(x1)-1;
j = 2:length(y1)-1;
k = 2:length(z1)-1;
V(i+1,j,k) = V(i-1,j,k);
V(i,j+1,k) = V(i,j-1,k);
V(i,j,k+1) = V(i,j,k-1);