@vv Your calculation shows great insight.  There is magic in the definition of H.  I was looking for something like that but was unable to derive it.  I have to admit that even after seeing it I still don't quite understand why it works, but that's OK since it gives me something to think about.

@Axel Vogt I don't see where you get the 0..10 range in z1.  That surely can't be correct.  Did you calculate that by hand or did you use a Maple library or something?

This is not a solution but only a sketch of a strategy which after a lot of effort will reduce the problem to a triviality.

The germ of the idea is observing that the integrand is constant on any hyperplane x₁ + x₂ + … + x₁₀ = r where r is a constant.  Therefore, changing the coordinates so that one of them points in the direction of that hyperplane's normal, the ten-fold integral reduces to an integral of a single variable in that direction.  The price we pay is having to deal with  geometric complications that arise from that change of coordinates.

I will demonstrate the method by solving the case of the triple integral.  The ten-fold integral is left as an exercise for the interested reader ;-)

Thus, let us look at the problem

where g is any function of one variable.  Change from the (x,y,z) variables to a new (u,v,w) variables defined through

x = u,
y = v - u,
z = w - v.

We see that x + y + z = w, therefore the integrand g(x+y+z) equals g(w), and the problem reducdes to that of the integration of a function of a single variable.

For the purpose of explanation, it helps to look at the more general problem

where f(x,y,z) is not necessarily of the form g(x+y+z).  To perform the change of variables, we examine the intersection of the cube [0,1]x[0,1]x[0,1] with the plane x+y+z=w and we see that there are three distinct types of intersections, depending on the value w.  When w is between 0 and 1, the intersection is an equilaterial triangle.  When w is between 1 and 2, the intersection is a hexagon.  When w is between 2 and 3, the intersection is again an equilateral triangle.  Thus, we slice the cube with the plane x+y+z=w as w varies from 0 to 3, integrate over each cross-section, and add up the results, which yields the vallue of the desired integral.

After a bit of geometric scrutinizing, we see that i₀ = i₁ + i₂ + i₃ + i₄, where

In partucular, if f(x,y,z) = g(x+y+z), each of these reduces to an integral of single variable g(w).  As an example:

which is approximately equal to 0.1589542911.

The attached worksheet has some extra details.

funky-integral.mw

@mskalsi dsubs is indeed the right tool for this.  I had forgotten about its existence.  Thanks for your note.

@mskalsi Ah, now I see what it is that you are asking.  I don't have a neat solution for you, but perhaps the attached worksheet will do.

zero-derivatives.mw

@Markiyan Hirnyk Yes, of course.  It is within the context of the solution of that system of PDEs that diff(u,x)=0.  Outside of that context u is just an undefined symbol.  That's why I prefaced my original comment by saying "I assume that you are loading the PDEtools package because you intend to solve a PDE".  I thought that my intention was clear.

@Markiyan Hirnyk A PDE is defined over a domain—that would be an open subset of the (x,y,t) space in this case.  One specifies initial and/or boundary conditions on the domain's boundary to arrive at a well-posed problem.

@Markiyan Hirnyk By "throughout" I meant "throughout the domain".

@Carl Love Your suggestion is a good one for a conventional eigenvalue problem, that is, finding ω so that det(A - ω I) = 0, where I is the identity matrix, or the more general det(A - ω C) = 0, where C is any matrix. The problem posted in the first worksheet in this thread, that is, eign.mw,  is neither. It is a nonlinear eigenvalue problem which seeks the zeros of det(A - C(ω)) = 0, where C(ω) is a nonlinear function of ω.

The problem in this thread's the second worksheet, that is, 1.mw, is of the form det(A - ω^2 I) = 0, which is amenable to your suggested method, as vv has observed.

@Carl Love Your suggestion is a good one for a conventional eigenvalue problem, that is, finding ω so that det(A - ω I) = 0, where I is the identity matrix, or the more general det(A - ω C) = 0, where C is any matrix. The problem posted in this thread is neither. It is a nonlinear eigenvalue problem which seeks the zeros of det(A - C(ω)) = 0, where C(ω) is a nonlinear function of ω.

This may be operating system dependent.  I tried your code on Maple 2015.2 on 64-bit Linux.  The memory usage stayed at 20.18 MB throughout the computation with n=4, 5, and 6.  I did not even "restart"  between the tests.

@iman I know that you don't work in astronomy.  My point was to make you think about the meanings of the numbers that you are working with.  The error in your calculations is not related to Maple at all.  The error already exists when you arrive at your matrix q.  As soon as you see numbers such as 10^17 and 10^(-118) mixed together in a matrix, you should stop, think, and ask yourself what it is that you are doing wrong and how to fix it.  Continue with calculations only after you have resolved that issue.

Think about this:  The q[1,1] entry in the matrix matrix q in your eign.mw worksheet contains the term 1.717631735 x 10^17.  Why is this number specified to only 10 digits of precision?  It's because you don't know the next digit -- the one that comes after "5".

What you don't know can be important.  If you insert a "1" just after the "5", the value of your number changes by ten million because it gets multiplied by 10^17.  So your number is accurate within plus or minus 10 million at best.  In view of this, what is the significance of the term 3.864000000^(-105)*omega^26 which also appears in q[1,1]?

To compute successfully, you should readjust the mathematical model to avoid a mixture of ridiculously large and ridiculously small numbers such as 10^17 and 10^(-108).  Don't expect to get far without doing that absolultely essentiual initial work.

@Sana Gull

Thomas Richard's code works just fine on Maple 11 and produces a result identical to that of the latest Maple.  I expect that it should work on your Maple 13 as well.  Try again.

@MatthewsonT I still don't see what it is that you are asking. As I wrote before, since you say you can do the problem on paper, it may be a good idea to show us exactly what it is that you do, providing all the details.  Then someone may be able to tell you how to do that same thing in Maple.

In the harmonic motion generated by the equation x'' + w^2 x = 0, the frequency is independent of the amplitude, therefore I don't quite understand what you mean by depicting one against the other.  Perhaps there is a hidden assumption which you have omitted.

Since you say you can do the problem on paper, it may be a good idea to show us what it is that you do.