@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

Of course.
The mapping (x,y) -> `if`(x+y<1, [x,y],[1-x,1-y])   applies the square [0,1]^2 onto a standard triangle similar to the mentioned Dirichlet distribution.

@Markiyan Hirnyk 

My code does not use incidence verification. It constructs effectively area preserving mappings just like Mathematica.

@Markiyan Hirnyk 

@Heitor 

It was just a joke, probably not a good one.

@torabi 

A general advise (at least for the future): it is a good idea to present mathematically the problem.  Otherwise what alternative method are you hoping for? It would imply a "reverse engineering" of your computations and probably nobody will want to spend his time for this. E.g. probably Runge-Kutta makes sense for the original problem, but not for the system you have presented.

@ecterrab 

What is the performance impact of Assume? After Physics:-Assume(...), Maple seems to use more memory.
It would be interesting to know why this operating mode was not included in the standard assume command (maybe via a kernelopts(assume) switch).