close all;

n = 4;
N=n*n;
Z = zeros(n);
I = eye(n);
I1 = diag([0 -ones(1, n-2) 0]);

T =   diag([1 4*ones(1, n-2) 1])
    + diag([-ones(1, n-2) 0 1])
    + diag([0 -ones(1, n-2) 1]);


Z_row = Z;
for i=1:n-2
    Z_row = [Z_row Z];
end
H(1:n,:) = [I Z_row];
H(N-3:N,:) = [Z_row I];

for i=2*n:n:n*(n-1)
    low_band = i-2*n+1;
    high_band = i+n;
    H(i-(n-1):i, low_band:high_band) = [I1 T I1];
end

b = zeros(N,1); b(15,1) = 10;

x = gmres(H,b);
r = reshape(x, [n,n]);
figure;
contourf(r,23);