@Carl Love 

convert/`if`  seems to be not documened.
What type of expression does it obtain?
P:=piecewise(x<7,x+1,x);
PP:=convert(P,`if`);
lprint(PP);   #     Error, invalid expression for eval

Edit. OK, you answered the question while I was asking, thank you  :-),

@abcd 

A random variable X (with values in T subset R^2) is considered uniformly distributed in the sense
Prob(X in D) = Area(D)/Area(T)
for each measurable subset D of T.

My sample generated for the ellipse satisfies this, while yours - for the unit disc (obtained from uniform and triangular distributions) does not.
You can verify this by counting the number of points in a disc (with radius 1/2 e.g.) not centered at the origin.
[You did not apply the definition, but a kind of statistical intuition].

@Markiyan Hirnyk 

Nevermind, it was pure curiosity.

@Markiyan Hirnyk 

Yes, but what about the Mathematica timing?

@Markiyan Hirnyk 

We do not know how Mathematica measures the timing.
So, I thing that the best way to compare would be to execute in a loop the procedure (maybe with distind input data) and compute the difference of cpu or real time (or even use a physical clock).
Also, the  graphical part should be excluded if we want to compare only the algorithms.

I also think that Mathematica will win only if the random generator is better, because the algorithm is too simple and taylored exactly for the problem.

@Carl Love 

It runs indeed very fast.
I wonder if a (compilable and) compiled version would increase the speed.
But if the bottleneck is the random generator ... that's it.

@Parham2016 

The boundary/initial conditions are needed for the uniqueness of the ODE solution. You can't modify them.
A similar situation would be the definition of sin via an ODE:
f'' + f = 0,  f(0)=0, f'(0)=1.

Changing to f(0)=1, f'(0)=0 ==> cos.

@Carl Love

The problem with such efficient procedures is that the logic (i.e. algorithm)
is sometimes obfuscated and the user has to rewrite (for him) some/many
lines in order to understand.
Here the situation is simple, except maybe
<<[1$n]> | A>.< <T[1]> | <T[2]-T[1]> | <T[3]-T[1]>>^+
but in other cases the procedure could be almost cryptic.

So, the ideal situation would be to have both versions posted.

@zia9206314 

Being entire, HeunT has only infinity as singularity, except when it reduces to a polynomial (and this happens for certain values of the parameters, see the help page).

I don't know if/when HeunT can be expressed via hypergeom.

@Markiyan Hirnyk 

So you have approximated 1/2 by 12/23.
Still, in the unit disc should be about 1333 points and there are 1448.
If you want to use this repartition instead of an exact one, it's OK.

@Markiyan Hirnyk 

You use improperly "unbased words" instead of "I need an explanation for ...".

In this case, f: [u,v]  |-->  [sqrt(u)*cos(v),  sqrt(u)*sin(v)]  has the Jacobian 1/2, so
area(f(D)) = 1/2 * area(D).

@Markiyan Hirnyk 

Because the mapping  [u,v]  |-->  [sqrt(u)*cos(v),  sqrt(u)*sin(v)]  is area preserving (modulo a constant),
which is not the case with  [u,v]  |-->  [u*cos(v),  u*sin(v)] .

PS. Thank you for the link. BTW, Robert Israel is one of the most remarkable scientists I ever met (unfortunately not personally); it's a pitty that he is not active any more in this forum.

@Markiyan Hirnyk 

You forgot to execute the first part which defines r.

You said first

... codes includes the verification of the incidence to the region (the tringle or ellipse).

wrongly and without any explanation.

@Markiyan Hirnyk 

And I understand that you have objections for the ellipse too. Please explain.

@Markiyan Hirnyk