%
%  Plot the function.
%
  r = 3;
  x = -r : 0.01 : +r;
  y = -r : 0.01 : +r;

  [ X, Y ] = meshgrid ( x, y );

  Z = ( X .* X + Y .* Y <= r * r ) .* ( X .* X );

  surf ( X, Y, Z, 'FaceColor', 'Interp', 'EdgeColor', 'Interp' )
  xlabel ( 'X' );
  ylabel ( 'Y' );
  zlabel ( 'Z=F(X,Y)' )
  title ( 'Z = X^2' )
%
  exact = 63.617251235193313079;
  area = pi * r * r;

  n = 1;

  for logn = 0 : 7

    x = - r + 2 * r * rand ( n, 1 );
    y = - r + 2 * r * rand ( n, 1 );

    i = find ( x .* x + y .* y <= r * r );
    n2 = length ( i );

    if ( 0 < n2 )
      result = area * sum (  x(i) .* x(i) ) / n2;
    else
      result = 0;
    end

    fprintf ( 1, '  %8d  %f  %e\n', n2, result, abs ( result - exact ) );

    n = n * 10;

  end
 
