plot3d(-abs(x-y) , x=-1..1, y=-1..1
    #, gridstyle=triangular
    , axes=boxed
 );

You can increase Digits or use more numpoints as well, 
see ?plot3d/option

What you do is to write Euler's integral as 2F1, but it is not
continuous in its parameters (else it would be 1/2 in that case)
and Maple does not do it.

To prove it: 

  eq3; combine(%);
  op(1,%);
  limit(%,k=infinity);
                                  0

So the integrand is 0 in the limit. Now justify interchanging
limit and integral to get your desired result.

Not that complete as the package as described in the linked paper above (especially they do have various methods), but the following fits my basic needs (and I am too 'lazy' for a full solution, that's on Maple, if they want to compete ...):

Download 102_computing_Gauss_integration_rules.mws (30 kb)

It has nothing to do with the package, perhaps you want the following code snippet:

  endless:=proc() while true do end do end proc; # endless waiting loop

  for i from 1 to 3 do
    try
      timelimit( 0.5, endless());
    catch "time expired":
      print(i, "can not wait"); #print(lastexception);
    finally
      print("going ahead ...", exp(i));
    end try;
  ``; # prints a blank line 
  end do;

                          1, "can not wait"
                      "going ahead ...", exp(1)

                          2, "can not wait"
                      "going ahead ...", exp(2)

                          3, "can not wait"
                      "going ahead ...", exp(3)

2F1 is just the more common notation for Maple's hypergeom with 2 parameters in the first and 1 in the second place.

A short sketch how to pass from series in the open unit disc to functions beyond it via integrals is given in

Abramowitz and Stegun

on p. 558 ff, online at www.math.sfu.ca/~cbm/aands/toc.htm

If you want to step deeper into those functions you may try

Lebedev

the beautiful book "Special Functions" comes as paperback for almost no money,
www.amazon.com/Special-Functions-Their-Applications-Lebedev/dp/0486606244
(though it is heavy stuff and real Math)

It is a bit hidden, You can do that for example with the package 'OrthogonalSeries':

Define your polynom from Taylor series:

  taylor(exp(-x)*sin(x), x, 11): 
  myP:=convert(%,polynom);

Convert it to the Legendre base and check it:

  OrthogonalSeries:-ChangeBasis(myP,LegendreP(v, x)):
  OrthogonalSeries:-ConvertToSum(%):
  convert(%,polynom):
  'myP=(%)';
  is(simplify(%));

                                 true

The same for Tschebyscheff:

  OrthogonalSeries:-ChangeBasis(myP,ChebyshevT(n,x)):
  OrthogonalSeries:-ConvertToSum(%):
  convert(%,polynom):
  'myP=(%)';
  is(simplify(%));
                                 true

If it is homework, then give credit to mapleprimes ...

clearer ... may be, I am too tiered :-) If it is not too ugly just upload your Excel sheet (you should 'zip' it to get it accepted), may be those days others or I find time to look into it

I googled a longer time to find out what you may mean. Besides some Matlab code for plotting from Excel the most reasonable I got was public.deltares.nl/display/HYMOS/05+IDF+Relations and may be that helps to state a more precise question. I guess it is to estimate a continous density from occurencies, where some specific model is taken - it seems, that in Hydrology this is the Gumbel distribution (or log-Weibull as a dictionary says). May be that also helps for formulating ... may be it is about estimating parameters (which may depend on time or time frames).

Played with it ... (but remember only vaguely, that Int(f(x), x=...)
constructs a procedure and I had cases, where this needs long - while
acer's approach is fine).

  Digits:=14;

  acer:="0.49597789486775255522822143063830103993131048383835*1e-10";

  convert(expr,rational): subs(699=a,%): h:=subs(700=a+1,%);
  select(has,h,a): g:=unapply(%,x);
  h/g(x): op(1,%)+op(2,%): f1:=unapply(%,x);
  h/g(x) - f1(x):  
  c:=denom(%);
  numer(%%): completesquare(%,x): f0:=unapply(%,x);

  'h=(f1(x) + f0(x)/c)*g(x)'; is(%);

                               true

This still does not have a reasonable slope, but the following does it,
only the last int has to be done numerically - using double precision:

  f:= x -> (f0(x)-f0(0))*g(x)/(a+1);
  F:= x -> eval((f0(x)-f0(0))*g(x),a=699);

  'int(f1(x)*g(x),x=0..1) + int(f0(0)*g(x),x=0..1)/c + (a+1)*Int(f(x),x=0..1)/c';
  'eval(%,a=699)';
  evalf(%);
                                           -10
                        0.49597789483778 10
  parse(acer);
                                           -10
                        0.49597789486775 10
                                    ^^^^

Even brute quadratures like (a+1)*student[simpson](f(x),x=0..1,20) do
result in the same figure, they agree only at the first 10 decimals
with the reference value 'acer'.

So try with higher precision:

  d:=50;
  Digits:=d;
  'int(f1(x)*g(x),x=0..1) + (int(f0(0)*g(x),x=0..1) + 
    Int(F,0..1, epsilon=0.5*10^(-d+2),digits=d+2))/c';
  eval(%,a=699):
  evalf[d](%);
  parse(acer);
  Digits:=14;

                                                             -10
      0.49597789483778839837287771760761560977019494673212 10

                                                             -10
      0.49597789486775255522822143063830103993131048383835 10


The different result occurs as well. Hm ...

you divide by 0 in theta = Pi/2

For me the task is a bit questionable, as it is at the limits
for usual floating point precision (say 14 Digits):

The integrand writes as f(x)*g(x).

g(x) = (a+1) * (2729/2728864913*x+1999899/1999901)^a, a = 699,
hence 700 * some_line^a, some_line ~ 1.

But (1-small)^699 can not be done reasonable in double precision.

f(x) is almost a constant of magnitude 'numerical_epsilon'.

Hence I am not astonished that the NAG method will not work.

And I doubt the input is 'exact' and using high precision can
not heal it. May be a reasonable scaling is desirable before
asking for the integral.


Some brute, but not that bad way, is to approximate f by its
Taylor polynomial (degree 2 or 3 should do, in x=1/2 if you
want some averaging), since f is almost constant.

Though it is always good to elaborate and understand answers of the
seniors like Robert Israel you may wish to re-think your way and goal
(i.e. I agree with Thomas Unger):

do the same for -x^2 on -1 ... 1 and moreover restrict the view to
the values by a command like 

  plot(-x^2,x=-1..1, y= 0 .. -eps)

with eps = 1e-5 or 1e-6 (may depend on Digits, version and interface)

This is not quite readable and it may be better, if you upload the question as a Maple file (try to use Classical interface, extension *.mws), which contains your formulae.

We use that for home office, it even blocks Excel books in the
new xml format (which is the standard for Ecxel 2007).

However it blocks only the possiblity for downloading those
attachments (giving an info about that).

Just tried for a test sheet: it can be sent. Opening this mail
shows the attachment, but downloading through mouse klick does
only offer 'attachment.ashx' which seems to be almost empty
(it is 1 KB, can see nothing in it).

Note that this may depend on your browser (mine is Firefox, it
behaves different from IE for OWA and one needs lower security
setting for IE to work there).

As this OWA is only a kind of interface to my actual mailbox
I can access it as POP3 (using Thunderbird) and there it was
no problem at all to download, detach and run the test sheet.

  int(int(2*x*sin(y),y=0..1),x=0..1);
#                             ^

             1 - cos(1)