Можете помочь с переводом от матлаба в с++ уравнение Пуассона?

% MATRPASS.M
%
% fi=matrpass(x,y,f,aH,bH,cH,aV,bV,cV,flag)
%
% Осуществляет решение двумерного уравнения Пуассона в декартовых (x,y)
% или циллиндрических (r,z) координатах с граничными условиями первого, второго или
% третьего рода (первый порядок аппроксимации) методом матричной прогонки
% с автоматическим выбором направления прогонки
%
% x,y - координатные векторы, задающие разностную сетку
% f - матрица размерности (length(x),length(y)), описывающая правую частиь
% (с обратным знаком)
% aH,bH,cH - матрицы размерность (2,length(x)), задающие коэффициенты a, b и c
% граничных условий на горизонтальных границах: a(x)U+b(x)dU/dy=c(x)
% (первая строка соответствует нижней границе, вторая - верхней)
% aV,bV,cV - то же, для вертикальных границ
% (первая строка - левая, вторая - правая)
%
% flag - переключатель системы координат:
% flag=0 - декартовая, flag<>0 - цилиндрическая
%
% [Cамарский А. А., Гулин А. В. "Численные методы математической физики",
% Научный Мир, Москва, 2000, с. 187-191]
%
function fi=matrpass(x,y,f,aH,bH,cH,aV,bV,cV,flag)
% размерность задачи
N1=length(x)-1;
N2=length(y)-1;
% Xpass - флаг прогонки по x, Ypass - флаг прогонки по y
if N1>N2
Xpass=0;Ypass=1;
else
Xpass=1;Ypass=0;
end;
% граничные условия разностной задачи
% (определяются из соответсвующих условий исходной дифференциальной задачи)
% на вертикальных границах
for i2=1:N2+1
% нижняя граница
denom=(aV(1,i2)*(x(2)-x(1))+bV(1,i2));
Hi1(i2)=bV(1,i2)/denom;
Mu1(i2)=cV(1,i2)*(x(2)-x(1))/denom;
% верхняя граница
denom=(aV(2,i2)*(x(N1+1)-x(N1))+bV(2,i2));
Hi3(i2)=bV(2,i2)/denom;
Mu3(i2)=cV(2,i2)*(x(N1+1)-x(N1))/denom;
end;
% на горизонтальных границах
for i1=1:N1+1
% левая граница
denom=(aH(1,i1)*(y(2)-y(1))+bH(1,i1));
Hi4(i1)=bH(1,i1)/denom;
Mu4(i1)=cH(1,i1)*(y(2)-y(1))/denom;
% правая граница
denom=(aH(2,i1)*(y(N2+1)-y(N2))+bH(2,i1));
Hi2(i1)=bH(2,i1)/denom;
Mu2(i1)=cH(2,i1)*(y(N2+1)-y(N2))/denom;
end;
% знаменатель второй разностной производной
% в направлении x (~квадрат шага)
for i1=2:N1
h1(i1-1)=(x(i1+1)-x(i1))*(x(i1)-x(i1-1));
end;
% ...в направлении y
for i2=2:N2
h2(i2-1)=(y(i2+1)-y(i2))*(y(i2)-y(i2-1));
end;
if flag
% знаменатель первой разностной производной
% в направлении x (для цилиндрических координат)
for i1=2:N1
ha(i1-1)=(x(i1+1)-x(i1-1))*x(i1);
end;
end

% задание матриц для прогонки по x
if Xpass
% вектора правой части системы матричных уравнений
F=zeros(N2-1,N1-1);
for i1=1:N1-1
F(1,i1)=(f(i1+1,2)+Mu4(i1+1)/h2(1))*h1(i1);
for i2=2:N2-2
F(i2,i1)=f(i1+1,i2+1)*h1(i1);
end;
F(N2-1,i1)=(f(i1+1,N2)+Mu2(i1+1)/h2(N2-1))*h1(i1);
end;
% коэффициенты матрицы оператора второй разностной производной по x
% необходимые в случае цилиндрических координат
if flag
for i1=2:N1
a1(i1-1)=1-flag*h1(i1-1)/ha(i1-1);
a2(i1-1)=1+flag*h1(i1-1)/ha(i1-1);
end
end
% матрица опрератора второй разностной производной по y
L=zeros(N2-1);
L(1,2)=1/h2(1);
for i2=2:N2-2
L(i2,i2-1)=1/h2(i2);
L(i2,i2)=-2/h2(i2);
L(i2,i2+1)=1/h2(i2);
end
L(N2-1,N2-2)=1/h2(N2-1);
% матрицы (диагональные) коэффициентов граничных условий
% левая граница
X1=diag(Hi1(2:N2));
% правая граница
X2=diag(Hi3(2:N2));
% векторы свободных членов граничных условий
% левая граница
M1=Mu1(2:N2)';
% правая граница
M2=Mu3(2:N2)';
end;

if Ypass
% вектора правой части системы матричных уравнений
F=zeros(N1-1,N2-1);
for i2=1:N2-1
if flag
F(1,i2)=(f(2,i2+1)+Mu1(i2+1)*(1/h1(1)-1/ha(1)))*h2(i2);
else
F(1,i2)=(f(2,i2+1)+Mu1(i2+1)/h1(1))*h2(i2);
end;
for i1=2:N1-2
F(i1,i2)=f(i1+1,i2+1)*h2(i2);
end;
if flag
F(N1-1,i2)=(f(N1,i2+1)+Mu3(i2+1)*(1/h1(N1-1)+1/ha(N1-1)))*h2(i2);
else
F(N1-1,i2)=(f(N1,i2+1)+Mu3(i2+1)/h1(N1-1))*h2(i2);
end
end;
% ма