@ebrahimina 

If Maple  does not use parallel algorithms for these, it will be very tough for you to implement them.

@ebrahimina 

I don't think that Maple can speed up significantly the computations for this precision.
You will need special tools and hardware. For example there are processors having the floating point arithmetic on 128 bits (==>  34 decimal digits) but it is not widely supported.

P.S. Maybe a better algorithm exists for your problem.

@s144117 

You must be aware that many functions (including hypergeom) are defined (in Maple and other CAS)  beyond the convergence domain (series or integral) using analytic continuation.
A very simple example:
f := hypergeom([1], [], x);
is defined by the series Sum(x^n, n=0..infinity) which converges for |x|<1. But its analytic continuation is f = 1/(1-x) defined for x in C \ {1}.
simplify(eval(f, x=7));
     - 1/6

@s144117 

1. It's not the integral; the series itself diverges (so, the integral is out of the question).

2. It's the ratio test: a series Sum(a[n], n=1..infinity) diverges if Limit(a[n+1]/a[n], n=infinity)  > 1.

restart;
Pn:=proc(n::posint, lambda, mu::Not({0,0.0}), P0, t::name)
local k,P;
P[0]:=P0;
P[1]:=(diff(P[0],t)+lambda*P[0])/mu;
for k to n-1 do  
    P[k+1] := (diff(P[k],t)+(lambda+mu)*P[k]-lambda*P[k-1])/mu
od;
P[n]
end:

Pn(5, 1/4, 1/2, exp(-t), t);
                            - 3 * exp(-t) / 32                               

@Zeineb 

Of course, but only after you define P[0](t), n, lambda, mu.

@sand15 

The OP mentioned "uniform distribution", so probably Prob(X=x) = 1/(k+1).
But anyway he does not seem to be interested in the answer.

I have a technical question, knowing that you are an expert in this kind of things.
1. I have never used thismodule. Trying the example in the help page, it crushed Maple 2018. In Maple 2017 it works. Do you know why? Bug?

2. You have used in this post the construction thismodule['aname']  to access the local variable aname using its global name.
Is this documented?

@student_md 

So, you want to define a sequence of matrices.
It would be better to use as index the position of the matrix in the sequence. I.e. C(k) is the k-th matrix (having the order 2^(k-1) in your case.

You need to define clearly H(k), c(k). With your notation for c_a = c_(m/2) it is not clear what is e.g.  c_3  in the sequence c_0, c_1,...,c_(m/2-1). That's because with your notation, c(m) is defined only for m = "a power of 2".

In other words, use consistent math notations. It will be easy then to obtain the maple code.

@Kitonum 

Of course, but the workaround often works when the locus is not known.
And after all your example is a (degenerate) ellipse.

@Kitonum 

But once we are aware of the (inherent!) problem, workarounds are possible. E.g.
plots:-implicitplot(f(x,y)=d+1e-5, x=0..6, y=0..6, gridrefine=5);

@student_md 

The main advantage is that you can edit it (or even create it with a text editor).

An .m file (internal format) loads faster, but your file is small and does not matter.

@peter2108 
You already have writebytes for this.

@ecterrab 

I just want to note that in the first integral (J), the change of variables was not done at all. J was simply rewritten using intat.

@Rouben Rostamian  

The EPS export for 3D is (and was) totally wrong (at least in Windows 64 bit).
plot3d(x^2+y^2, x=-1..1,y=-1..1);
(I have posted twice the result in the past.)