@Carl Love 

I don't understand. You said (correctly) "The problem is not caused by using a variable whose name begins with underscore" and then you provide a workaround exactly for such variables (only).

restart;
assume(Z1::integer);
sys:={x <> (  (1/2)*Pi+Z1*Pi ), x < infinity, -infinity < x};
expand(thaw(solve(subsindets(sys, suffixed(_), freeze@``) ,{x}))); #error

Also, the freezed variable is without assumptions; so, a simpler woraround would be to simply remove the assumptions which are anyway ignored now.

@Carl Love 

Yes, you are right; and for other types too, e.g. positive.

@Preben Alsholm

@digerdiga 

Using such strange names for variables (for an outsider) your code looks as being obfuscated and very hard to follow.
Why don't you isolate the problem and use normal (short) variables?

@Annonymouse 

subsindets(P, specindex(integer,x), _XX)

replaces each x[n]   having n::integer  by _XX(x[n]);
Then,
remove(has, %, _XX);
removes the polynomials containing _XX

@nm 

It seems that Maple ignores the singular points. You can simply exclude them from Carl's solution:

sing := remove(has,[singular(f)],infinity);

@Adam Ledger 

The function is discontinuous; you may add option discont.

Choosing these coords=...  appears pointless, unless you have a good reason for it.

AFAIK there is no user-defined storage functions for Arrays.  A possibility is to use storage=sparse and write the index function; it would be interesting to see the efficiency in this case.

@Markiyan Hirnyk 

It would be nice to tell us your ideas. E.g.
1. A multi-dimensional Van Aubel theorem.
2. A complex program for animated flexible polyhedra.
3. A project to generalize Gaifullin's results.
4. Something else?
 

@mmcdara 

0<x<1, 0<y<1 is not necessary because we are going to integrate the product of densities.

E.g. for a=4/3 Maple gives;

So, it's not unreasonable to expect from PDF (when a=4/3) the expression

piecewise(z<=0, 0, F(z))
but it seems that you are ok with the present status.

@mmcdara 

I expect from PDF to return at least the (double) integral (not necessarily a closed form, of course, even if such a closed form exists). This can be easily done for the toy example because Maple is able to solve Z < z:

SolveTools:-SemiAlgebraic({x^2 + y^2 + a*x*y<z},[x,y], 'parameters'=[a,z]);

(not to mention the case when a is fixed, e.g. a=1/2 or a=3).

So, I don't agree about "Not really a severe limitation".

@Alger 

It is not clear what you are trying to do. If you want to find an orthogonal system in a function space, the simplest way is to start with a known one e.g. cos(nx), sin(nx), in L²([0,2*Pi], 1)  and use the Fact above for some h.

@Rouben Rostamian  

The fact applied here is simple:

If u_n(x) are orthogonal in L²([x1,x2], w(x))  (here w is the weight) and h : [t1, t2] --> [x1,x2] is a diffeomorphism
then   u_n(h(t))  are orthogonal in L²([t1,t2], w(h(t))*|h'(t)|).

@minhthien2016 

Thank you for the new recurrence.
Have you tried to replace the new coefficients in the provided answer?

@Mariusz Iwaniuk 

On my systems (Maple 2017 or 2018, Win 7, 64 bit) it works in all situations: 1D/2D, Document/Worksheet.