function value = quad_fun ( n  )

  fprintf ( 1, 'Compute limits\n' );
  spmd
    a = ( labindex - 1 ) / numlabs;
    b =   labindex       / numlabs;
    fprintf ( 1, '  Lab %d works on [%f,%f].\n', labindex, a, b );
  end

  fprintf ( 1, 'Each lab estimates part of the integral.\n' );

  spmd
    x = linspace ( a, b, n );
    fx = f ( x );
    quad_part = ( b - a ) * ( fx(1) + 2 * sum ( fx(2:n-1) ) + fx(n) ) ...
      / 2.0 / ( n - 1 );
    fprintf ( 1, '  Approx %f\n', quad_part );
  end

  quad = sum ( quad_part{:} );
  fprintf ( 1, '  Approximation = %f\n', quad )

  return
end
  
