if i call your set of unknowns by unk, then there are 22, yes.
but indets(eqns,atomic): indets(%): % minus unk; nops(%); shows
{s1, s2, v0, s0, vQX, vQY, v1, v2, RT1_FR, RT2_FR, RL1_LS, RL2_LS, PaQv}
so your input produces 13 new variables ... what is s1, s2 etc?
May be a careful input may solve your problem?

That would be fine ... but they would have done, if they think it time
to release it.
And it might not help you - why 50 digits (I think Bailey has special
software for quadrature of high precision)? Plotting shows the problem
is at the boundary.
Some thoughts ... as you do not state all you know for your special
case and I see
- 2 variables (could there be more? - that may be a show stopper)
- coeffs sum to 0 in (x,y)=(0,0), so re-scaling possible?
- there is only 1 mixed term
- staring at the powers: one can reduce to degree 2
Int(Int(sqrt(R)/xi^(99/100)/eta^(1410/1411),xi = 0 .. bX^100),
eta = 0 .. bY^1411)/176375000 where bX, bY are the bounds and
R:=1562500+7582500*eta+9199089*xi^2+7582500*xi+9199089*eta^2
-33563178*xi*eta
Then the problems for the integrator are at 0. Note that order of
integration numerical could make a difference, if the powers are
of different size (as in your case, bounds^m --> 0 with m large)
and one of the bounds is not very close to 1.
Now i looked in Gradsteyn et al ... 2.260.1 has a recursive formula
for Int( x^m * srqt(R), x). One might try to write the denominators
as Taylor series (may be after xi -> 1 - mu etc) and try to write
down a recursion, so the problem becomes an infinite series in
applying the outer integration numerical (plus some stopping, I do
not have a good feeling for convergence rates ...).
But for higher dimension this approach will be quite hopeless and
it reminds me a bit to calculate a (bi-variate) cumulative normal
where correlations ~ +-1 cause problems, but can be settled in the
2 dim case (only) by recursive computations.
But it's just an idea, not a solution.

you just type the formulae in a sheet ...
and for all software new to you: contact the manual for first steps ...
sounds sarcastic, I know, but either i do not understand your question ... or you have not looked into the manual or the only help (just mark in with the mouse and hit F1 key

could you say what you want to do exactly?

There is the possibility to tell Maple to use a default type for the variables.
The other - quite lame - thing I sometimes do is to work with an editor or just
Win Word (which often is easier than to look up commands, especially if one wants
arrays in readable form or other formatings): copy your definitions to the editor
and replace "t" by "float t" using the editor's ability.
But the proper way is to look at some of the examples at the online help.

This is not a function of P_m, since that is your integration variable. If the others are constant that is a real number (or undefined) and I doubt that Maple plots it - what is your input for that?

Even if I accept licence checking as a good thing that is either
an unacceptable behaviour or a bonelazy answer you got. Or both.
Of course I have files which are created in future ... I do not
use my PC for Maple alone.

It is not quite clear for me which constants you use and it may be you cut off
to early (where the exp term is still contributing), but I would do it this way:
Your integral writes as Int(erfc(b*sqrt(xi/(xi+1))) *exp(-a*xi)*a,xi=0..infinity)
for b = B^2, a = 1/gamav.
Here xi/(xi+1) approaches 1 from below for large xi and beyond xi = 1/epsilon-1
differs at most by epsilon from 1. Then split the integration at this point into
Int(a*erfc(b*(xi/(xi+1))^(1/2))*exp(-a*xi),xi = 0 .. 1/epsilon-1) +
Int(a*erfc(b)*exp(-a*xi),xi = 1/epsilon-1 .. infinity);
The last is erfc(b)*exp(-a*(1/epsilon-1)) and the first should be done by working
numerical.
May be you refine that suggestion.

since nobody answers: you may add your hardware and operating system
and may be graphical interface (and personally i just would reinstall
instead of doing anything else ... and de-install first ... and if
you do it try to avoid the blank in the directory name, overwrite it
by Maple10 or similar within the installation dialog ... and if you
have and old OS like WinME - do not use the standard interface)

would you mind to post that example? evalf(Int(f,t=0..1)) uses the NAG library i think ...

randomize():
B:=[stats[random,normald[0,1]](N+1)]:
W:=[stats[random,normald[0,1]](N+1)]:
Z:= zip((x,y) -> evalf( rho*x+sqrt(1-rho^2)*y ),B,W): # or evalhf?

Then Z and B will do, where rho is the desired correlation.
In Maple 10 it is a bit different, but much (!) faster.

Once i translated Ooura's double exponential integration into
Maple (find it uploaded, but i have neither cleaned up that
port nor tested it carefully). That method is quite robust.
So you can use Maple's new ability to generate a DLL and use
that as a function directly in Maple. This might give you
some of the speed improvements you are looking for.
The additional advantage is: you can test values against Maple
if you have some doubts.
For distributions usually a Gauss-Legendre method also works
well (using fixed numbers of sub-divisions), 16 point suffice.
What i can not see is your chi^2, i thought it contains an exp
in the integrand? Just a last thought: when i try to fit some
pdf against data (to get parameters) i pass to the logarithm
(and refine later if neccessary). If you want to fit against
a cdf then i can imagine troubles (i am not sure whether this
is a stable approach, especially if you are not absolutely
convinced that your cdf is the correct one).
PS: as i do not recognize how to attach files use that link
http://www.axelvogt.de/axalom/maple/Ooura_DoubleExp_Integration.mws

the only thing i can imagine is that Maple methods can be
worked out ... having not traced through the procs it is
clear that they involve external callings for NAG routines
what exactly are you looking for?