% 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