Assuming you saw something:
>> K=4*diag(ones(4,1))-diag(ones(1,3),1)-diag(ones(1,3),-1) K = 4 -1 0 0 -1 4 -1 0 0 -1 4 -1 0 0 -1 4 A1=[K,-eye(4,4),zeros(4,4),zeros(4,4)] A2=[-eye(4,4),K,-eye(4,4),zeros(4,4)] A3=[zeros(4,4),-eye(4,4),K,-eye(4,4)] A4=[zeros(4,4),zeros(4,4),-eye(4,4),K] A=[A1;A2;A3;A4]; A(1:8,1:8) >> A(1:8,1:8) ans = 4 -1 0 0 -1 0 0 0 -1 4 -1 0 0 -1 0 0 0 -1 4 -1 0 0 -1 0 0 0 -1 4 0 0 0 -1 -1 0 0 0 4 -1 0 0 0 -1 0 0 -1 4 -1 0 0 0 -1 0 0 -1 4 -1 0 0 0 -1 0 0 -1 4 >> b=[130,20 10 80 60 0 0 60 40 0 0 40 20 0 0 20]'
>> tvec=A\b tvec = 52.0212 33.5727 28.1318 36.9939 44.5121 34.1379 31.9606 39.8439 31.8894 26.5061 25.7288 30.4212 16.5394 14.2682 14.0273 16.1121
tmat=reshape(tvec,4,4) surfc(tmat); %colorbar % rotate, see contours etc.
Continue to Laplace and Poisson equations