抛物型差分（二维-1）

// ImpMx.m
function ImpMx = getImpMx()
clear;clc;
 
M = 4;
N = 3;
Up = zeros(1,M);
Down = Up;
Lf = zeros(1,N);
Rt = Lf;
rx = 0.1;
ry = 0.3;
r = rx + ry;
ImpMatrix = diag( ones(1, M * N)) * ( 1 + 2 * r);
size( ImpMatrix)
b = zeros(M * N, 1);
for k = 1:N
    for j = 1:M
        rowIndex = index2pos(j ,k ,M);
        if( j - 1 == 0)
            b( rowIndex) = b(rowIndex) + Lf(k);
        else
            ImpMatrix( rowIndex, index2pos(j -1,k, M)) = -rx;
        end
        if ( j + 1 == M + 1)
            b(rowIndex) = b(rowIndex) + Rt(k);
        else
            ImpMatrix( rowIndex, index2pos( j+1,k,M)) = -rx;
        end
        if( k - 1 == 0)
            b( rowIndex) = b(rowIndex) + Down(j);
        else
            ImpMatrix( rowIndex, index2pos(j,k -1,M)) = -ry;
        end
        if( k + 1 == N + 1)
            b( rowIndex) = b(rowIndex) + Up(j);
        else
            ImpMatrix( rowIndex, index2pos( j,k+1,M)) = -ry;
        end
    end
end
 
 
ImpMx = ImpMatrix;
 
//index2pos.m
function pos = index2pos( m, n, N)
pos = m + (n - 1) * N;

注： 1.未给出具体的解，但构造了与一维类似的矩阵，其一行最多有5个非零元素
    2.所得矩阵非3对角矩阵，但也有解法，故考虑其他的差分格式
    3.ADI格式下次再给出，即交替差分格式