function x = diophantine_5d_nonnegative ( a, b )

%*****************************************************************************80
%
%% diophantine_5d_nonnegative finds nonnegative diophantine solutions.
%
%  Discussion:
%
%     We are given a Diophantine equation 
%
%       a1 x1 + a2 x2 + a3 x3 + a4 x4 + a5 * x5 = b
%
%     for which the coefficients a are positive integers, and
%     the right hand side b is a nonnegative integer.
%
%     We are seeking all nonnegative integer solutions x.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 May 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer a(5): the coefficients of the Diophantine equation.
%
%    integer b: the right hand side.
%
%  Output:
%
%    integer x(5,k): solutions to the equation.
%
  if ( length ( a ) ~= 5 )
    error ( 'diophantine_5d_nonnegative: a does not have length 5.' )
  end

  if ( any ( a <= 0 ) )
    error ( 'diophantine_5d_nonnegative: a has a nonpositive entry.' )
  end

  if ( ~ i4vec_is_integer ( a ) )
    error ( 'diophantine_5d_nonnegative: some a is not an integer.' );
  end

  if ( sum ( a ) <= 0 )
    error ( 'diophantine_5d_nonnegative: a does not have a positive entry.' );
  end

  if ( b < 0 )
    error ( 'diophantine_5d_nonnegative: b is negative.' );
  end

  if ( ~ i4_is_integer ( b ) )
    error ( 'diophantine_5d_nonnegative: b is not an integer.' );
  end

  d = i4vec_gcd ( a );
  if ( mod ( b, d ) ~= 0 )
    error ( 'diophantine_5d_nonnegative: b is not divisible by gcd ( a ).' );
  end

  x = [];
  y = zeros ( 5, 1 );

  a5 = [ 1; a(5) ];
  b5 = b;
  [ d5, v5, w5, cmin5, cmax5 ] = diophantine_2d_nonnegative ( a5, b5 );
  for c5 = cmin5 : cmax5
    y(5) = v5(2) + c5 * w5(2);

    a4 = [ 1; a(4) ];
    b4 = b5 - a(5) * y(5);;
    [ d4, v4, w4, cmin4, cmax4 ] = diophantine_2d_nonnegative ( a4, b4 );
    for c4 = cmin4 : cmax4
      y(4) = v4(2) + c4 * w4(2);

      a3 = [ 1; a(3) ];
      b3 = b4 - a(4) * y(4);
      [ d3, v3, w3, cmin3, cmax3 ] = diophantine_2d_nonnegative ( a3, b3 );
      for c3 = cmin3 : cmax3
        y(3) = v3(2) + c3 * w3(2);

        a2 = [ 1; a(2) ];
        b2 = b3 - a(3) * y(3);
        [ d2, v2, w2, cmin2, cmax2 ] = diophantine_2d_nonnegative ( a2, b2 );
        for c2 = cmin2 : cmax2
          y(2) = v2(2) + c2 * w2(2);
          a1 = a(1);
          b1 = b2 - a(2) * y(2);
          if ( mod ( b1, a1 ) == 0 )
            y(1) = b1 / a1;
            x = [ x, y ];
          end
        end
      end
    end
  end

  return
end

