function [ a, f ] = boundary ( nx, ny, node_num, node_xy, indx, ib, nunk, a, f )

%% BOUNDARY modifies the linear system for boundary conditions.
%
%  Modified:
%
%    17 May 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NX, NY, controls the number of elements along the
%    X and Y directions.  The number of elements will be
%    2 * ( NX - 1 ) * ( NY - 1 ).
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, real NODE_XY(2,NODE_NUM), the coordinates of nodes.
%
%    Input, integer INDX(NODE_NUM), gives the index of the unknown quantity
%    associated with the given node.
%
%    Input, integer IB, the half-bandwidth of the matrix.
%
%    Input, integer NUNK, the number of unknowns.
%
%    Input, real A(3*IB+1,NUNK), the NUNK by NUNK
%    coefficient matrix, stored in a compressed format.
%
%    Input, real F(NUNK), the right hand side.
%
%    Output, real A(3*IB+1,NUNK), has been adjusted for boundary conditions.
%
%    Output, real F(NUNK), has been adjusted for boundary conditions.
%

%
%  Consider each node.
%
  node = 0;

  for row = 1 : 2 * ny - 1

    for col = 1 : 2 * nx - 1

      node = node + 1;

      if ( row == 1 | row == 2 * ny - 1 | col == 1 | col == 2 * nx - 1 ) 

        i = indx(node);
        x = node_xy(1,node);
        y = node_xy(2,node);
        [ u, dudx, dudy ] = exact ( x, y );

        jlo = max ( i - ib, 1 );
        jhi = min ( i + ib, nunk );

        for j = jlo : jhi
          a(i-j+2*ib+1,j) = 0.0;
        end

        a(i-i+2*ib+1,i) = 1.0;

        f(i) = u;

      end

    end
  end