function x = rref_solve ( a, b )

%*****************************************************************************80
%
%% rref_solve() uses reduced row echelon form (RREF) to solve a system.
%
%  Discussion:
%
%    A linear system A*x=b is given, where
%    A is an MxN1 matrix, possibly singular.
%    b is an Mx1 right hand side, or MxN2 collection of right hand sides.
%    x is the desired N1x1 solution, or N1xN2 collection of solutions.
%
%    The right hand sides are assumed to be consistent, that is,
%    each right hand side is assumed to be a linear combination of 
%    columns of A.  If this is not so, the procedure will still produce
%    a result x, but it will not satisfy the equation A*x=b. 
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    28 September 2018
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A(M,N1), a matrix.
%
%    real B(M,N2), one or more right hand sides consistent with A.
%
%  Output:
%
%    real X(N1,N2), one or more solution vectors.
%
  [ m1, n1 ] = size ( a );

  [ m2, n2 ] = size ( b );

  if ( m1 ~= m2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'RREF_SOLVE - Fatal error!\n' );
    fprintf ( 1, '  B does not have the same number of rows as A.\n' );
    error ( 'RREF_SOLVE - Fatal error!' );
  end

  ai = [ a, b ];

  [ airref, k ] = rref ( ai );

  x(1:n1,1:n2) = airref(1:n1,n1+1:n1+n2);

  return
end

