@ThU: If I do what you suggest, and then subsequently inserts a new section, then the problem seems to reappear.

@Preben Alsholm: There is indeed no reason to expand all sections.

@sunit: Happy to hear that you got it working.

@sunit: I am afraid that that did not clarify matters further. I think my answer is still as before. Sorry.

@sunit: In order to avoid using subs, you could perhaps define a new procedure:

myDiff := (f) -> convert(diff(f,t),diff):

Applying it to lambda(t - tau_1), gives the usual result:

myDiff(lambda(t - tau_1));

Applying it to your more complicated expression gives:

newFunction :=
   +lambda[1](t-tau_1, alpha*t-alpha*tau_1, alpha^2*t-alpha^2*tau_1)
   +lambda[2](t-tau_1, alpha*t-alpha*tau_1, alpha^2*t-alpha^2*tau_1)*alpha
   +lambda[3](t-tau_1, alpha*t-alpha*tau_1, alpha^2*t-alpha^2*tau_1)*alpha^2:
myDiff(newFunction);

Without being sure at all, I guess the convert part of myDiff does not work here because the arguments in lambda1, etc., are not independent. For comparison, consider the following two study cases:

convert(diff(f(x - a,y    ),x),diff);
convert(diff(f(x - a,x - b),x),diff);

Perhaps it would be profitable to express lambda1, etc., as functions of only two arguments: the independent ones t - tau_1 and alpha, for the three arguments as presently given may be written as t - tau_1, alpha*(t - tau_1), and alpha^2*(t - tau_1).

What is the purpose of the quantities A,B,C,E,F,G? I cannot see them figuring in PolynomialSystem below.

@sunit: No need to apologize.

Does tau_1 mean t1? Are lambda1, etc., functions of three arguments?

Also, could you provide the example you give in a form where I can copy-paste it to Maple?; at present it is an image.

@Carl and John_Dirichlet: I see.

Just being curious: Why not use factor to put it on the following much more nice form?

factor(w^2/(w^4+2*w^2+1));

@Markiyan Hirnyk: Thanks for correcting my mistake. It would seem that I have to be more precise. Unfortunately, by training I am not a mathematician, but just a physicist, with all the mathematical sloppiness, I guess, that often accompanies that. I would really love to be able to take some formal course in distribution theory, for my "math problem", as you correctly put it earlier, remains frustratingly unresolved.

Concerning the wiki-page: Just underneath Eq. (2) there is the following warning (which I first noted after posting above): "However, despite widespread use in engineering contexts, (2) should be manipulated with care, since the product of distributions can only be defined under quite narrow circumstances."

Update concerning the Wikipedia Eq. (2): Others believe it to be entirely correct.

@Markiyan Hirnyk: Thanks again for your feedback. And thanks for having a go with the CPV calculation.

Perhaps we are talking past each other on the issue of Fourier transformation. Let the issue be, if it is ok with you.

This may become embarrasing for me (if I am fundamentally mistaken, that is), but I have always thought that the n-dimensional Dirac delta function was just the product of n one-dimensional Dirac delta functions of the n different arguments, see for instance Eq. (2) here, at least if these arguments are all independent. That was why I asked as I did.

@Preben Alsholm: Thanks for correcting my mistake. Somehow I had a hunch that my answer was not completely correct. But just to be sure: the problem with the posters code originates in the code line I claim, right?

@Markiyan Hirnyk: Thanks for your interesting input.

I am fully aware that the Fourier transform I am looking for is a generalized one, i.e., a distributional one in the sense of generalized functions. This follows from the simple fact that Heaviside(x) + Heaviside(-x) = 1, so that the Fourier transformation of Heaviside(s^2) + Heaviside(-s^2), where s^2 = t^2 - x^2 - y^2 - z^2, will give (2*pi)^4 times the 4-dimensional Dirac delta function.

But the result Mathematica delievers surprises me. Just to be sure: it is 4*pi^2 times the 4-dimensional Dirac delta function, right? If indeed so, then it seems completely wrong to me, for two reasons:
1.) Only a constant has Fourier transform proportional to the Dirac delta function. Heaviside(s^2) does not satisfy that condition because it vanishes outside the light cone in Minkowski space.
2.) My pen and paper calculations yield a result in the form a*DiracDelta(E,px,py,pz) + b*P(1/k^4), where a and b are some constants, where k^2 = E^2 - px^2 - py^2 - pz^2, and where P denotes the Cauchy principal value.

Anyway, what I was (as well) asking for was whether Maple has similar higher-dimensional calculation features in versions later than 17. Looking at your solution using VectorCalculus:-int, DotProduct, etc., I guess not, or what?

@Carl Love: I agree, using op is kind of a 'dirty' solution. Thumbs up to your nice solution, again one from which I learn new about Maple.

@HelenS: Happy to hear that.