function obj = fe_system( rho, pts )
    % The first argument, 'rho', appears in Poisson's
    %           __ 2
    % equation: \/  u = rho
    %
    % The second argument, 'pts', is a 3 x N array
    
    if size( pts, 1 ) ~= 3
        error( 'expecting a 3 x n matrix of points but got a %d x %d matrix', ...
            size( pts, 1 ), size( pts, 2 ) );
    end
    
    n = size( pts, 2 );
    from_alias = zeros( n, 1 );
    to_alias = zeros( n, 1 );
    ni = 0;

    % For any point 'pts(i)' which is set to -Inf, we record
    % the value in the tables to_alias and from_alias.
    
    % The variable 'ni' stores the number of entries to be found.
    
    for i = 1:n
        if pts(3,i) == -Inf
            ni = ni + 1;
            from_alias(ni) =  i;
            to_alias(i) = ni;
        end
    end
 
    % We will create a system Mu = b and solve for u once
    % an equation has been defined with each point
    M = zeros( ni, ni );
    b = zeros( ni,  1 );
    % Initially, we have updated none of the interior points
    %  - we can update each interior point once
    updated = zeros( ni, 1 );

    obj = struct( 'rho', rho, 'n', n, 'ni', ni,         ...
        'points', pts, 'updated', updated,              ...
        'from_alias', from_alias, 'to_alias', to_alias, ...
        'M', M, 'b', b );
end

