@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.

@Carl Love 

FAIL/FAIL and even FAIL/0 or 0/FAIL returns FAIL, so the previous version seems to be ok (with the typo corrected).

If d=0 why have you decided to return 1? It should be FAIL (or remain 0/0, hence error) because in this case the definition of the conditional probability is different, so not handled by the formula.

@tomleslie 

It was just a remark. I agree with you of course that the results of the two parsers should be the same.

Using complex numbers the proof is even shorter:

P[5]:=P[1]:
for i to 4 do Q[i] := P[i] + (P[i+1]-P[i])*exp(-I*Pi/4)/sqrt(2) od:
is( Q[3]-Q[1] = I* (Q[2]-Q[4]) );

        true

n dots (n>=2) surrounded by spaces (which is probably the preferable syntax) are interpreted as  ..  in both 1D and 2D.
seq(-3 ... 3);  seq(-3 .... 3);

@Markiyan Hirnyk 

The typo is obvious and the line was shown.

Conditioned:= proc(
   event::{relation, set(relation)},
   cond::{relation, set(relation)}, 
   rv::
      {name= RandomVariable &under Statistics:-RandomVariable,
       {set,list}(name= RandomVariable &under Statistics:-RandomVariable)
      }
) 
uses St= Statistics;
local 
   Event:= `if`(not event::set, {event}, event),
   #Cond:= `if`(not cond::set, {cond}, event),   #typo
   Cond:= `if`(not cond::set, {cond}, cond), 
   RV:= `if`(not rv::{set,list}, [rv], [rv[]]) 
;
   if nargs=3 then 
      (Event,Cond):= 
         subs(lhs~(RV)=~ St:-RandomVariable~(rhs~(RV)), [Event, Cond])[]
   fi;
   #union used rather than intersection because Probability considers a set
   #of relations the represent the intersection of the corresponding events. 
   St:-Probability(Event union Cond)/St:-Probability(Cond)
end proc:

Conditioned({2*X+4*Y <= 1}, {4*X+2*Y <= 1}, {X = Uniform(0, 1), Y = Uniform(0, 1)});

      FAIL

@Markiyan Hirnyk 

After correcting the mentioned typo it will return FAIL.
P.S. If Probability works for not necessarily independent vars, then Conditioned would be ok.