@abcd 

You are right. Both our approaches are correct.
The reasons I suspected that you are wrong are:
1. I am a 'pure' mathematician. I am not used to see  r + dr  etc. For me it's easier to compute a Jacobian
(see a previous intervention).
2. The Triangular distribution was misleading. I would have used the (equivalent in this case) Power distribution.
3. I had the visual impression that the density is higher off the center.

All the best, and sorry for my mistake.

What are  f^n(x_0) and df^n(x_0)/dx ?

@Carl Love 

Probably a good idea (from an esthetical point of view) would be
`convert/If`:= ...

convert(... , 'If');

PS. Is such functional style considerably more efficient?

@Carl Love 

piecewise(x<1,10,x<2,20,30);

is equivalent to

`if`(x<1,10,`if`(x<2,20,30));

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