And what is not so simple here?

@Markiyan Hirnyk

@one man

It would be nice to say something about the method used in the worksheet.
What "universal parametrization" have you used? Otherwise, the reader must use a "reverse engineering"  technique to find out.

If the surface is already parametrized, one may use something like this:

CP:=(v,w)->[v[2]*w[3]-v[3]*w[2], v[3]*w[1]-v[1]*w[3], v[1]*w[2]-v[2]*w[1]]: #Cross Product
IP:=(v,w)->v[1]*w[1]+v[2]*w[2]+v[3]*w[3]: # Inner Product
UV:=v -> v/~sqrt( IP(v,v)); #Unit Vector

r:=[s*cos(t), s*sin(t), s^2]:          # Example, paraboloid (parametric)
Nr:=UV( CP( diff(r,s), diff(r,t) ) ): #Unit Normal vector to r
Er:=simplify(  r + 1/2*Nr) ;           # Equidistant to r


p1:=plot3d( r, s=1 .. 2, t=0..2*Pi ):  p2:=plot3d( Er, s=1 .. 2, t=0..2*Pi ):
plots[display](p1,p2);






Edit. OP's animation using this method:

The problem is that Maple gives the result without assumptions, even if for a complex k, the integral may not exist (converge) e.g. for k = I.

@m3rcur1al 

You may use the Explore command and play with the parameters, as in the example:

Explore(fsolve(c1*exp(x)-c2*sin(x), x), parameters=[c1=1..15, c2=-3..13]);

You may prefer to replace fsolve with a plot.

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