Matlab fehlerhafte Systemmatrix für die Randbedingung Neumann?

Hallo Leute,
ich habe folgendes Problem. Ich muss eine elliptische DGL in einem beliebigen rechteckigem Gebiet mit Neumann Randbedingungen in Matlab programmieren und lösen.
Nur mit meiner programmierten Systemmatrix funktioniert das nicht und ich bekomm unendlich große Werte. Habe das Problem auch schon zum Teil herausbekommen (glaube ich).

Hier erstmal der Programmabschnitt für die Matrix:

h=1/3; % Festlegung der Schrittweite h   
k=1/4; % Festlegung der Schrittweite k
hk=k/h; % k faktoriell in h abgebildet
%---Anfang Systemmatrix Neumann---%
%---Anfang Blockdiagramm---%
d1=diag((2+2/hk^2)*ones(n,1));
r1=-diag(1*ones(n-1,1),1);
l1=-diag(1*ones(n-1,1),-1);
d1=d1+r1+l1;


d1(1,2)=-2; 
d1(n,n-1)=-2; 
%---Ende Blockdiagramm---%


%---Anfang Seitendiagonale rechts---%
r2=-diag((1/hk^2)*ones(nm-n,1),n);


for kk=1:n
    r2(0+kk,n+kk)=-1/hk^2*2;
end
%---Ende Seitendiagonale rechts---%


%---Anfang Seitendiagonale links---%
l2=-diag((1/hk^2)*ones(nm-n,1),-n);


for kk=1:n
    l2(nm+1-kk,nm-n+1-kk)=-1/hk^2*2;
end
%---Ende Seitendiagonale links---%


mat=d1;


for kk=1:m-1
    mat=blkdiag(mat,d1);
end


mat1=mat+r2+l2;


clear mat;
%---Ende Systemmatrix Neumann---%
disp(mat1)
mat1=sparse(mat1);
%;

Die Matrix, die erstellt wird, sieht genauso aus, wie die, die wir in den Vorlesungen auf dem Papier erstellt haben.
Wenn ich die zwei "-2" im Block einen anderen Wert gebe, bekomme ich vernünftige Zahlen.

5.5556 -2.0000 0 0
-1.0000 5.5556 -1.0000 0
0 -1.0000 5.5556 -1.0000
0 0 -2.0000 5.5556

Ich hoffe, es ist nur ein kleiner Fehler, den mir ein gübter Matlab'er/erin beantworten kann.
Erstmal vielen Dank, ich werde morgen wieder vorbeischauen.