Roughly you want

6.6e18906*n*((-2.302585093+ln((-6.*n+7.*n*m+4.*n*t+8.-10.*t-10.*m)/(-2.*n+2.*n*m+n*t-3.*m-4.*t+3.)))/(.5e6990-.1e6990*n)
 ^(7447/10000)/m)^(10000/2753)*m +
690.7755279-.3e3*ln((-6.*n+7.*n*m+4.*n*t+8.-10.*t-10.*m)/(-2.*n+2.*n*m+n*t-3.*m-4.*t+3.)) = 0.

As it was said 6.6e18906 is likely to kill numerics - where you would have to supply concrete
values for m=taum and t=taut to work at all.

And symbolically I have doubts there are usable solutions.

Best would be to think about appropriate scaling of variables, which you have to do before
arriving at your final task.

Just as thoughts

u[x](x)*(y/delta(x)-2*y^2/delta(x)^2+y^3/delta(x)^3+
   B*H*(1-3*y^2/delta(x)^2+2*y^3/delta(x)^3)/delta(x))/
  (3*B*H/delta(x)+2);
 
 Int(%^2, y=0 .. delta(x)): value(%):
 
 simplify(%);

                2                   2                         2  2
         u[x](x)  delta(x) (delta(x)  + 11 B H delta(x) + 39 B  H )
   1/105 ----------------------------------------------------------
                                               2
                           (3 B H + 2 delta(x))

find the pdfs attached maple10_sheets_pd.zip

You may try to upload them here, as zipped files

# use 'exact' values
[m = 0.211, c = 0.80, r = 0.338, e = 1];  params:=convert(%, rational);

eq:=H(z)*(1-m*H(0)^2*(1+z)^3/H(z)^2)-2*e*sqrt(1-c^2)*sqrt(r)*H(0) = 0;

For H(0) = 0 you only have H = 0, otherwise:

eq*H(z): collect(%, H(z));
[solve(%, H(z))];
eval(%, params); 
sols:=simplify(%);

which gives 2 solution for H(z), H(0) not zero.

Due to your eq you finally have to check, that H(z) <> 0,
map(solve, sols, z); gives the chance for z = -1.

eval(eq, z=-1);
eval(%, params);

relates H(0) and H(-1), where I stop, lacking further ideas.


Edited:

For the given parameters I get H(z) = 0, I think (H(0) Real). 
May be, that you need the parameters more precisely.

The double integral in the special case corr=0=mean, variance=1 is a product of
normals (as already said) and you curve p = G(y1)*G(y2) can be parametrized as
y2 = H( p/G(y1) ), G = cdf normal, H = G^(-1) the inverse

G := t -> 1/2+1/2*erf(1/2*t*2^(1/2)), H := y -> RootOf(-erf(_Z)-1+2*y)*2^(1/2),

Q:= (y,p) -> H(p/G(y)); # gives the curve, like a hyperbola

To be defined the input has to be between 0 and 1 and this defines 'bounds', as
they are gives by the equation p = ..., since G is between 0 and 1 and p = G(y)
gives boundary := -1.64485362695147 (use Digits = 15), beyond it becomes Complex.

You will see, that these are the 'asymptotics' in plotting it:

  'Q(y1, 0.05)';
  plot([%, boundary - y1, 0*boundary], y1 = boundary .. 6, boundary .. 6,
    labels=["y1","y2"]);

Of course y1+y2 is not constant on it, since Q is like a hyperbola.

Now you want the minimum of y1 + y2 = y1 + Q(y1,p) (justify: the *same* y1 !)

Plotting it suggests y1 ~ -1.0, as you said, so set up diff(t+Q(t, 0.05), t);
fsolve(%, t=-1); giving -0.760068575155508. And it turns out that y2=Q(...)
is the same (i.e. a fixed point ... strange), so the minimum is ~ -1.520

That is also formally correct:

If the pair is on the curve, then it must satisfy p = G(t)^2, t:=y1=y2 and it
should be a minimum, for which one substitutes that for p:

  'eval(diff(t+Q(t, p), t), p=N(t)^2)';
  simplify(%) = '%';

                       /d               \|
                   0 = |-- (t + Q(t, p))||
                       \dt              /|        2
                                         |p = N(t)

So it is an extremal point (I leave it to you to show that D@@2 is positive).

Hence one has: a minimum for y1+y2 is y1=y2, such that N(y1)^2 = p.

There is no global maximum (since y1+Q(y1) ~ y1 + const) and a local maximimum
at the left side in y1 = boundary.



PS: Perhaps best is to look into books / literatur on portfolio optimization and
'value at risk'. Certainly the working with normal distributions is off reality.
And using higher dimensional integrals may be not the best idea (even in dim=2
you will not get nice numerics if correlation is not zero, but beyond dim=3 that
will be expensive). However I never cared enough for that field, so I can not
recommend something besides what you would find anyway through Google or looking
at financial forums.

I am not very used to Statistics and after some times I forget definitions, always.

But if I get it right you want p = Int( Mpdf(x), x in cube ) where the range is a
cube, which starts in the 3 infinities and ends in corners x1,x2,x3.

But that is like asking  Function(x) - p = 0, a level or zero set of a function,
which for 3 variables 'generally' is a surface.

And I am not sure, whether the package is designed for that - what do you expect
as output or result?


The other point is:

If that is a probability function to be integrated, then usually towards boundary
the integral becomes quite flat approaching 0 or 1.

Numerically (if asking for 1 point) it is often much better to solve the equation
log( F(x) ) = log( small ).


But may be I got all wrong ...


PS:

for your last post at
http://www.mapleprimes.com/questions/129055-Quantile-And-Multivariate-PDF#comment129063
you still will have p = Function(xx[1],xx[2], x[3]), so the above remark on level
set applies, that is asking for a surface in R^3.

Likewise here one can do:

ShortHand:= my= J(s)/V(s);
Prepare:=   isolate(%, J(s));

eval((Ls+R+1/Cs)*J(s)=V(s), Prepare);
isolate(%, my);
eval(%, ShortHand);

                         J(s)         1
                         ---- = -------------
                         V(s)             1
                                Ls + R + ----
                                          Cs

That also works for L*s, C*s instead of Ls, Cs

(f@@3)(x); 
expand(%);

                              f(f(f(x)))

see the help about @

For the mentioned example

Hope it is usefull (though it should be cleaned up a bit)
and fast enough( about 1/2 second for that example,
if using Maple's slow BesselK).

Ooura_intde_oscillat.mws

Using Ooura's intdeo (for oscillating factors) and Maple's BesselK I get 
    -0.109408897877946e-3 
(or -0.109408897877997e-3 depending on treating the tail towards infinity)

while the reported Maple result through NAG + Ooura's BesselK0 is
    -0.109403613457678e-3

intdeo is his implementation for 'double exponential integration'
for oscillating factors to an integrand g(x) * sin(omega * x + theta)
where one provides the function g and frequency = omega = dz,
integrals are taken over a range ] constant ... infinity [

I think you need not care so much for memory (if you are in Maple),
the system should care. Have you tried 'forget'?

For BesselK 0: you can try evalhf'ing Ooura's code (or even compile it).

The usual way for such tasks is to approximate F =BK0(..)/BK0(...).
If your a and rho are constant that could be done once.

Since that F is like exp * 1/k one could try a ChebPade approximation
(I think, that is what Ooura does in some sense).

If one 'factors' it into partial fractions over the Reals, then Maple can
find Int( F*cos ) symbolically in terms of ?Si, ?Ci and that evaluation
will be fast (and would avoid special methods for high oscillations).

PS: I do not understand, why replying to a post sometimes works
and sometimes creates a new branch in the thread ... :-(
Can some admin shift it place, please ...

At least 1/2, but I can not prove it (but plotting and experiments 'show' it).

A 'crazy' task. 

First one thinks "if denominator not zero, one is ok" which would be Pi/3 ...

Then (working out taylor coeffs in 0 and using gfun) I used the solution for 

  {(2-54*t)*diff(y(t),t)+(4*t-27*t^2)*diff(y(t),`$`(t,2))-6*y(t), y(0) = 1, D(y)(0) = 3}

with y(z^3)*Pi = g(z), but that works only for z < 1/2.


Then I tried an analytic solution ( K = the integrand):

  K(z,t); expand(%): subs(cos(t) =c, %): simplify(%): 
  simplify(%, size): subs(c=cos(t), %);

Now Change(Integral, t = arccos(xi), xi) and completesquare gives 'value' the
chance to return a (complicated) result,

it is summing over the roots of a polynomial of degree 6 for a summand of
log*rational.

Inspecting shows, that in z=1/2 there _might_ be troubles, but they resolve :-(

And numerical testing give me the impression, that this solution is not correct
for z <= 1/3*2^(2/3) [that is, where the first way has imaginary solutions]

I give up ... 

Some observations, while waiting for more, specific informations about your task.

Let me write w for your variable dz, being positive and large.

Your kmin is positive and small, so you do not need 'piecewise' in the definition
for your function f.

Also I would avoid to involve 'evalhf' in its definition, the NAG routine will do
that for you, if possible.

Roughly your kmax turns out to be of magnitude ~ 30*Pi and your task writes as

  Int(F(k)*cos(k*w),k = kmin .. kmax) for F = k -> BK0(k*rho1)/BK0(k*a)

So it is a finite cosine transform.

If w is large, then this may be a highly oscillating integrand - already for the
cosinus term. And if F is oscillating as well (BK0 = BesselK ?) this may even be
worse. But that heavily depends on the concrete setting (if F has poles ...).

Changing coordinates w = t*Pi/k it becomes

  Pi/w * Int( F(t*Pi/w) * cos(t*Pi), t = 1/2 .. kmax/Pi*w )

so the upper bound is roughly 30*w and w is large.

There are special methods for these highly oscillating integrals, it may be that
NAG's method as 'working horse' is not appropriate.

TrapRule:=proc(f,N,a,b)
local k,sum1,h,i:
h:=(b-a)/N:
sum1:=0:
for i from 1 by 1 to N-1 do
#                    ^
sum1:=sum1+f(a+i*h):
end do:
k:=(h/2)*(f(a)+2*sum1+f(b)):
#                   ^
return (k):
end proc;

TrapRule(exp   ,4,-3,1); # use a function, not an expression
#           ^^^
           1/2 exp(-3) + exp(-2) + exp(-1) + 1 + 1/2 exp(1)

 evalf(%);
                             2.887249173