Point 1:
About your last sentence: "Final goal: select an arbitrary located circular domain for integration and a bivariate distribution with different variance in x and y"
There will be no close form expression of this integral.
To getone you must integrate over ellipses whose half axes are equal to the standard deviations of x and y.
In the more general case of non centered correlated components x and y, the integration has to be done done within elliptic domains centered at the mean x-y point, and whose axes are aligned along the eigenvectors of the variance matrix.

See wiki  item Cumulative distribution function

Point 2:
A main result in Probability Theory (see reference above) is that (I only consider here your simple case of centered uncorrelated components with equal standard deviations sigma) is that the probability that the random vector < X, Y > (whose PDF is your formula (1)) has a square norm less than r__mx² is related to the ChiSquare distribution with 2 degrees of freedom:

restart:
with(Statistics):
Z := RandomVariable(ChiSquare(2)):
Probability(Z <= (r__mx/sigma)^2) assuming r__mx > 0, sigma > 0
                             /        2 \
                             |   r__mx  |
                      1 - exp|- --------|
                             |         2|
                             \  2 sigma /

Point 3:
In case you woulf be interested in integrating a mulnormal PDF over ,o, canonical domains (ellipsoids), it's good to know there exist very specific fast algorithms to do this.

EDITED: I read @Carl Love's reply after I sent this reply. Maybe, @Carl Love,  you had in mind the same questions as me?


Assuming you have this expression

fcn1 := sum(N[i]*z[i]^2, i);
                        -----           
                         \              
                          )            2
                         /    N[i] z[i] 
                        -----           
                          i             

and want to compute by hand its derivative with respect to N[i],  what do you expect to get?

Let me be more explicit and assume fcn1 is written this way

fcn1 := J -> sum(N[i]*z[i]^2, i=1..J);
       J             
     -----           
      \              
       )            2
J ->  /    N[i] z[i] 
     -----           
     i = 1           

For, let's say, J=3 one gets

fcn1(3)
                       2            2            2
              N[1] z[1]  + N[2] z[2]  + N[3] z[3] 

and, for all value i in {1, 2, 3}, the derivative of fcn1(3) wrt N[i] is equal to z[i]^2.

So my question: "Would you expect Maple doing this:

                                      2 
Differentiate fcn1(J) wrt N[i]) = z[i]  

?"

If so, and by cheating a little, you could do this

Download Is_this_cheating.mw


I said before I was cheating because  I do not differentiate the sums but their generic terms, which happens to be the same thing because the expressions are indeed sums.

Had you define fcn1 as a product and not a sum, the result would have been wrong:

fcn4 := product(N[i], i):
Diff(fcn4, N[i]) = Derive(fcn4, N[i]);
                         /,--------'    \    
                         |   |  |       |    
                    d    |   |  |       |    
                  ------ |   |  |   N[i]| = 1
                   dN[i] |   |  |       |    
                         \    i         /    

(In the same logic as before the expected result is the product of all the N[k] where k takes any integer value but i.).
The_product_case.mw

Another way to understand I was cheating is that what  I did doesn't give the correct answer (assuming it is z[k]^2) to this question "What is derivative of fcn1 wrt N[k] where k is a positive integer?"


 

A few workarounds

b := add(rand(1 .. 6)(), i=1..3):
b := `+`(seq(rand(1 .. 6)(), k=1..3)):

Download Funny.mw

Using the Statistics package makes the things clearer

Download Statistical_equivalence.mw

eq1 := 2*A*B/(A*m+m)+C^2+D^2;

SET := A=1/alpha:
eval(eq1, SET):
map(simplify, %):
eval(%, isolate(SET, alpha));
                      
                           2 B       2    2
                      ----------- + C  + D 
                         /1     \          
                      m | -  + 1|          
                         \A     /          

One_Way.mw

The matrix M depends on several quantities:

• lambda
• m0
• a
• c
• N

One could write a procedure which build M from these  five quantities and computes its characteristic polynomial wrt some auxiliary quantity eta.
I didn't go in this direction and kept the things simple.
In the attached file you will find:

• a quite simple way to build M for arbitrary values of N,
• the computation of its characteristic polynomial  thanks to LinearAlgebra:-CharacteristicPolynomial,
• the computation of this same characteristic polynomial by explicitely computing the determinant ofmatrix
  M-eta*Identity(N+1) wrt the last row of this matrix.

CP.mw

Your notation FT{d(U(x, t)/dt}  = d(FT{U(t))} / dt, is confusing (see @Carl Love for another interpretation)

• if U(x, t) represents a function in the physical space, then its Fourier transform is:
  • some function H(xi, t) if you do a Fourier transformation wrt x,
  • some function H(x, tau) if you do a Fourier transformation wrt t.
• So do we read the rhs of your equation H(xi, t) or H(x, tau) ?
  Given you write d(FT{U(t))} / dt, I'm leaning to diff(H(xi, t), t) (otherwise you should have written diff(H(x, tau),tau) ), meaning the Fourier transform is likely taken wrt x: fourier : U(x, t) --> H(xi, t)

restart:
with(inttrans):

Notations:
A function named by a lower case letter applies in the physical space.
The same letter capitalized denotes the fourier fransform of a "physical" function.

I_set := fourier(u(x, t), x, xi) = U(xi, t);
               fourier(u(x, t), x, xi) = U(xi, t)

We_have := 'fourier'(diff(u(x, t), t), x, xi) = fourier(diff(u(x, t), t), x, xi);
          / d                \    d                         
   fourier|--- u(x, t), x, xi| = --- fourier(u(x, t), x, xi)
          \ dt               /    dt                        

By_notation := rhs(We_have) = eval(rhs(We_have), I_set);
            d                             d          
           --- fourier(u(x, t), x, xi) = --- U(xi, t)
            dt                            dt      
   
Then := 'fourier'(diff(u(x, t), t), x, xi) = rhs(By_notation)
                  / d                \    d          
           fourier|--- u(x, t), x, xi| = --- U(xi, t)
                  \ dt               /    dt         

Note that a Fourier transform wrt t gives

I_set := fourier(u(x, t), t, tau) = U(x, tau);
              fourier(u(x, t), t, tau) = U(x, tau)

'fourier(diff(u(x, t), t), t, tau)' = eval(fourier(diff(u(x, t), t), t, tau), I_set);
                / d                 \                  
         fourier|--- u(x, t), t, tau| = I tau U(x, tau)
                \ dt                /                  

... the rhs has nothing to do with diff(U(x, tau), tau) !

fourier.mw

restart:
kernelopts(version) 
      Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895 

a := 2*m:
b := 3*m:

# This can be done differently in more recent versions

'a*b' = ``(a)*``(b);
                       a b = (2 m) (3 m)

``(a)*``(b) = a*b
                                        2
                       (2 m) (3 m) = 6 m 

`a*b` := ``(a)*``(b)  = a*b
                                            2
                    a*b := (2 m) (3 m) = 6 m 

Solving numerically the Lane-Emden equation with theta__0(xi) being itself the numerical approximation of the Emden equation returns this error message

Error, (in dsolve/numeric/make_proc) ode system is singular at the initial point

Note that this has nothing to do with the fact the Emden solution is solved numerically:

EmdenE := unapply( rhs(dsolve({ode1(1), theta__0(0) = 1, (D(theta__0))(0) = 0}, theta__0(xi))), xi);

ode2E := n -> diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = -EmdenE(xi)^n*xi^2;
ic2  := gamma(0) = 0, D(gamma)(0) = 0;

LaneEmdenN := dsolve({ode2E(1), ic2}, numeric, parameters=[alpha])
Error, (in dsolve/numeric/make_proc) ode system is singular at the initial point

The reason why dsolving/numeric the Lane-Emden equation comes probably from the nullity at xi=0 of the denominator of its exact solution:

LaneEmdenE := dsolve({ode2E(1), ic2})
                (cos(xi alpha) xi alpha - sin(xi alpha)) _C2
    gamma(xi) = --------------------------------------------
                                     xi                     

         /    /     2    \   2\                       
         \2 + \alpha  - 1/ xi / sin(xi) - 2 xi cos(xi)
       + ---------------------------------------------
                               2            2         
                 xi (alpha - 1)  (alpha + 1)          

See at the end of this answer how to et the mw file.

So the idea to rewrite the initial conditions this way

ic2 := gamma(epsilon) = 0, D(gamma)(epsilon) = 0;

where epsilon is aimed to be a small positive quantity (please to the first point in  Remarks above).

The mw file will be attached as soon as this will be possible (the site experiences a technical problem presently).

Meanwhile here ifs the code (the important lines whowing how to use the numerical solution of the Emden equation are yellow highlighted)

restart
local gamma:

ode1 := n -> (diff(xi^2*(diff(theta__0(xi), xi)), xi))/xi^2 = -theta__0(xi)^n:
EmdenN := dsolve({ode1(1), theta__0(0) = 1, (D(theta__0))(0) = 0}, theta__0(xi), numeric):
plots:-odeplot(EmdenN, [xi, theta__0(xi)], xi=0..1):

f := proc(z)
  if not type(evalf(z),'numeric') then
    'procname'(z);
  else
    eval(theta__0(xi), EmdenN(z));
  end if;
end proc:
  
ode2 := n -> diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = -f(xi)^n*xi^2:
ic2  := gamma(epsilon) = 0, D(gamma)(epsilon) = 0:

LaneEmdenN := dsolve(eval({ode2(1), ic2}, epsilon=0), known=f, numeric, parameters=[alpha])
Error, (in dsolve/numeric/make_proc) ode system is singular at the initial point

# define an extra parameter epsilon
LaneEmdenN := dsolve({ode2(1), ic2}, known=f, numeric, parameters=[alpha, epsilon]):

# with epsilon=0
LaneEmdenN(parameters=[2.356, 0]):
plots:-odeplot(LaneEmdenN, [xi, gamma(xi)], xi=0..10):
Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

# with epsilon=1e-6 
LaneEmdenN(parameters=[2.356, 1e-6]):
plots:-odeplot(LaneEmdenN, [xi, gamma(xi)], xi=1e-6..10):  # works well

Remarks:

• In your file you used 2 conditions at xi=0 and another one at xi=xi_I: that's one to many.
• I used ode1 := n -> ... (the same for ode2) to give you the possibility to change in a simple way the value of n.
• I used the 'parameters' option for the same flexibility reasons.

Of course you can set n, alpha and epsilon to numerical values BEFORE dsolving the equations:

data := {n = 1, alpha = 2.356, epsilon=1e-6}:

ode1   := (diff(xi^2*(diff(theta__0(xi), xi)), xi))/xi^2 = -theta__0(xi)^n:
ic1    := theta__0(0) = 1, (D(theta__0))(0) = 0:
EmdenN := dsolve(eval({ode1,ic1}, data), numeric);

ode2 := diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = f(xi)^n*xi^2;
ic2  := gamma(epsilon) = 0, D(gamma)(epsilon) = 0;

f := proc(z)
  if not type(evalf(z),'numeric') then
    'procname'(z);
  else
    eval(theta__0(xi), EmdenN(z));
  end if;
end proc;

LaneEmdenN := dsolve(eval({ode2(1), ic2}, data), known=f, numeric)

Last but not least,
You will see in the mw file that the exact and numerical solutions of the Lane-Emden solutions are quite different.
I haven't been able to fix this in a satisfactory way.

To overcome the present technical difficulties Mapleprimes faces, you can download the Maple file this way:

• edit this link  DropBox,
• copy/paste it on your navigator,
• download the file (LaneEmdenNumeric.txt),
• change its name into LaneEmdenNumeric.mw,
• open LaneEmdenNumeric.mw with Maple.

Maybe this could help

For odd orders only:
The site is in French but the Maple code is in English :-)
http://villemin.gerard.free.fr/Wwwgvmm/CarreMag/aaaMaths/Tapis.htm

For odd orders and  few even orders
(Still in french)
Go here http://carres-magiques.com/partie-3-nouvelle
Then search Annexe 23

and click on "Cliquez ici pour obtenir les programmes Maple"
At last unzip the package

s := [{f0(t, r) = 0, n(t) = n(t)}, {f0(t, r) = 0, n(t) = 0}, {f0(t, r) = 0, n(t) = 0}]:

pattern := n(t)=n(t):
select(has, s, pattern)
                 [{f0(t, r) = 0, n(t) = n(t)}]

pattern := n(t)=0: 
select(has, s, pattern) 
[{f0(t, r) = 0, n(t) = 0}, {f0(t, r) = 0, n(t) = 0}] 

For this kind question randpoly is a useful tool.
But you can also write your own generator to create equations of a specific form.
For instance
 

possibilities.mw
 

Explanation of the parameter radfavor.
Replacing

t := `if`(rand(0..3)()<=2, -1, 1);

by

t := `if`(rand(0..19)()<=18, -1, 1);

will generate a radical in 95% of the cases.
If you write instead

t := `if`(rand(0..19)()=0, -1, 1);

a radical will be generated only 1 out of 20 times.

restart
B := Matrix(9$2, symbol=m):
Bd:= CodeTools:-Usage( LinearAlgebra:-Determinant(B) ):
memory used=0.82GiB, alloc change=0.60GiB, cpu time=17.64s, real time=11.98s, gc time=9.07s

# approximately the number of symbols used to represent Bd

L := length(Bd)
                            40642561
# printing Bd on a 400 charaters per line display will require
# this number of lines

iquo(L, 400) + `if`(irem(L, 400)=0, 0, 1)
                             101607

# And thus a pocket book (50 characters per line and 60 lines
# per page) of about this number of pages

round(evalf(L / 50 / 60))
                             13548

So how can you even read this expression?

On might think using the package ImageTools, providing the image has an appropriate format and is of good quality.
But it will be quite difficult to extract the contours of the points, take their centers, extract the correct scales (did I miss comething?)

My advice is to use a digitizer and do the things manually (some can be use online, for instance https://plotdigitizer.com/app)

... is to compute the determinant X of some symbolic matrix A of elements x_{i, j}, and next to replace in the expression of X each  x_{i, j} by the corresponding expression in B (x[i, j]=B[i, j])
The corresponding lines are in brown courier font.

Download determinant.mw
 

EDITED:
Following the (judicious) remark @acer  made, here is a simple say to get the determinant of B1 (not the one of Omeg):

detH  := Determinant(H):
detB1 := Determinant(B1):
detB  := detH * detB1:

Download determinant_edited.mw


As matrix Omeg := B + idn8 = H.B1 + idn8 (BTW idn8 := IdentityMatrix(8)) and because H:=Matrix(8, 8, symbol=h), (thus H is formally equivalent to the matrix X used in the "trick"), I think the "trick" can be applied to compute Determinant(Omeg).
Of course you must be aware of  @acer's remark if you have to evaluate Determinant(Omeg) for some instance of H.

UPDATED

What loop are you talking about? (shouldn't have been here)

Percentage of variation of peak heights
 

A_ref := eval(A, lambda = 1.3015):
x_ref := fsolve(diff(A_ref, x), x=-2..0);
y_ref := evalf(eval(A_ref, x=x_ref));

                         -0.7765307466
                          1.920953647

data_1 := [beta=0.1, Q=1.3015, lambda=0.9986];
B_1    := eval(B, data_1):
x_1    := fsolve(diff(B_1, x), x=-2..0);
y_1    := evalf(eval(B_1, x=x_1));

PercentageOfVariation := (y_1-y_ref)/y_ref*100

           [beta = 0.1, Q = 1.3015, lambda = 0.9986]
                         -0.8959052730
                          1.752467429
                          -8.770967392

data_2 := [beta=0.3, Q=1.3015, lambda=0.9986];
B_2    := eval(B, data_2):
dB_2   := diff(B_2, x):
x_2    := fsolve(diff(B_2, x), x=-2..0);
y_2    := evalf(eval(B_2, x=x_2));

PercentageOfVariation := (y_2-y_ref)/y_ref*100

           [beta = 0.3, Q = 1.3015, lambda = 0.9986]
                          -1.088722291
                          1.547856941
                          -19.42247313