% First, a simple case similar to that in Iskander's book (p. 307)
imax = 5; 
jmax = 9; 
bc = [1 1 0;1 2 0;1 3 0;1 4 0;1 5 0;1 6 0;1 7 0;1 8 0;1 9 100;...
    2 1 0;2 9 100;3 1 0;3 9 100;4 1 0;4 9 100;...
    5 1 0;5 2 0;5 3 0;5 4 0;5 5 0;5 6 0;5 7 0;5 8 0;5 9 100];
    
% Second, a finer grid case
imax = 9;
jmax = 17;
bc = [1 1 0;1 2 0;1 3 0;1 4 0;1 5 0;1 6 0;1 7 0;1 8 0;1 9 0;1 10 0;1 11 0;...
      1 12 0;1 13 0;1 14 0;1 15 0;1 16 0;1 17 100;...
    2 1 0;2 17 100;3 1 0;3 17 100;4 1 0;4 17 100;...
    5 1 0;5 17 100;6 1 0;6 17 100;7 1 0;7 17 100;8 1 0;8 17 100;...
    9 1 0;9 2 0;9 3 0;9 4 0;9 5 0;9 6 0;9 7 0;9 8 0;9 9 0;9 10 0;9 11 0;...
      9 12 0;9 13 0;9 14 0;9 15 0;9 16 0;9 17 100];

% Third, arbitary imax and jmax
imax = 5; jmax = 9;
nbc = 2*imax+2*jmax-4; 
[X Y] = meshgrid(1:jmax,1:imax); % prepare for boundary nodes numbering
mdl = (Y-1)*jmax + X;   % 1D indexing of boundary nodes
mdl(2:end-1,2:end-1) = 0;   % only boudary nodes have none zero indices
bc = zeros(nbc,2);  % 1D bc
bc(:,1) = mdl(mdl~=0); 
bc(end-imax+1:end,2) = 100;
K = imax*jmax; % total number of nodes

% build FD equations Ax = B.
A = eye(K)*-4;
for k = 1:K % for each row
    if k-1 > 0  % node at left arm
        A(k,k-1) = 1;
    end
    if k+1 < K  % node at right arm
        A(k,k+1) = 1;
    end
    if k-jmax > 0   % node at upper arm
        A(k,k-jmax) = 1;
    end
    if k+jmax < K   % node at lower arm
        A(k,k+jmax) = 1;
    end
end

B = zeros(K,1);

[I J] = size(bc);

if J == 3 % bc in the form of (i, j, v)
    bc = [(bc(:,1)-1)*jmax+bc(:,2) bc(:,3)]; % bc changed into (k, v)
end

for i = 1:I % boundary condition enforcement
    m = bc(i,1);
    A(m,:) = 0; A(m,m) = 1; B(m) = bc(i,2);
end

x = A\B; % solve for x
Phi = reshape(x,jmax,imax);
figure;image(Phi'); axis image; set(gcf,'color','w');
figure;contour(Phi',20); grid on; set(gcf,'color','w');