@Carl Love 

Of course it matters whether  x = x(beta,...)  because dx/dbeta etc will appear in the final result.

The problem is not correctly formulated, so an answer is not possible.
You must say which are the indepentent variables and which are constants.
Write for example r(beta,varphi) etc, if this is the case.

@Carl Love 

Do you know by chance what is the exact meaning of INTERVAL(a..b, c..d)  and INTERVAL(x, a..b)?
The help page for evalr mentions these but nothing more, so the user must guess a lot.

@tomleslie 

The ode here (OP's and my example) is in an extended sense, i.e. y(x) is continuous in [a,b] but the equation y'(x) = f(x, y(x)) is only for x in (a,b); so, y is not necessarily differentiable in a or b.

Note that in mathematics even more general situations are studied, e,g. y(x) is absolutely continuous, hence differentiable only a.e. (=almost everywhere, that is with the exception of a negligible set). See the wikipedia article on Carathéodory's existence theorem).

Consider the example
sqrt(x) * y'(x) = 1,  y(0) = 0.

The ode is of course for x>0.
Here the solution exists and is unique; it is continuous in [0, oo) and differentiable in (0, oo).

@mmcdara 

I get

f:=exp(-x^2-y^2/4-z^2/9):
a:=sqrt(k): b:=2*sqrt(k): c:=3*sqrt(k):
elli:=[x=a*r*sin(u)*cos(v), y=b*r*sin(u)*sin(v), z=c*r*cos(u)]:
simplify(int(eval(f,elli)*a*b*c*r^2*sin(u), u=0..Pi,v=0..2*Pi,r=0..1)) assuming k>0;

limit(%, k=infinity); # check

@radaar 

If someone is asking you:
I have a problem, I tried all the methods and they don't work. Can you help me?
What would be your answer?

@cduston 

I think that you miss the general picture. It is not easy to compute the residue if the order of the pole is not numeric.
In this case a numerical computation is not possible (actually the presence of any parameter makes a numerical integration impossible).
And we do not have an algorithm (the residue command simply uses (truncated)  series.
E.g. how would you compute  residue( exp(sin(z^2 - z^m)) / sin(z^n), z=0)  ? You will have to consider many possibilities for m and n.

(Not to mention that even in your example, n must be an integer, otherwise z=0 cannot be a pole).

@Rouben Rostamian  

The solution sol is actually a distribution. To obtain a classical one, f must be of course subject to some restrictions (otherwise e.g. the integral could diverge).

@Kaio 

Change is confused because you use the same symbol s in different contexts.
You can use

IntegrationTools:-Change(int((diff(f(s), s))^2, s = L .. ss), s = S*L, S);
    

and then
eval(%, ss=s);

@Kitonum 

I don't understand how you have obtained  `-`(a,b) = - a
It should be a-b, unless you have something like:

local `-`;
`-`:= (x,y) -> :-`-`(x);

@Magma 

1. Yes, all C[i,j] < 2^53. But if you are interested only in the positiveness of A^k this restriction can be relaxed.

2. If you want A^k to have positive entries, why do you need GF(2,n)? 
When I have mentioned Modular, I meant the use of a very large modulus m i.e. working in Zm. The role of m is similar here with the role of 2^53 for floats.
If I understand correctly your problem, the use of floats should be enough.

@Magma 

Search for randperm or type

?randperm

Or, simply type

combinat[randperm](5);

You will find it.

So, you just want some "cosmetic" rearrangements of the terms. Unfortunately (for some users) in Maple it's very difficult to do this; (and sometimes only using inert functions or Typesetting). Most users are happy to see the mathematical answers; the final aspect (e.g. for publication) is done by hand and is anyway very subjective.
But you already know the answer. Why not simply check it versus Maple's answer?  (Your answer uses different names for some constants).