I think you want something like this:

myCT := proc(L::{list, Vector})
   local locs, nL, n;   
   if type(L,list) then nL:= nops(L)   
   else nL:= LinearAlgebra:-Dimension(L)      
   end if;
   locs:= select(t -> type(L[t],posint), [$1..nL]);
   n:= nops(locs):
   if n = 0 then [[],[]]
   else 
      [[seq(L[locs[j]], j=1..n)],
       [locs[1],seq(locs[j]-locs[j-1],j=2..n)]]
   end if;
end proc;

Note that if you want to apply this to the output of Statistics:-Sample(Poisson(...),...)
you can leave that as a Vector but you do need to convert the entries to integers.  Thus:

L := map(round,Statistics:-Sample(Poisson(0.5), 10000));

time(myCT(L));

            0.031

Hints:

1) P(n,x) are orthogonal on the interval [-1,1] with int(P(n,x)^2, x=-1..1) = 2/(2*n+1)

(int doesn't seem to know this if you use LegendreP or orthopoly[P], but see the help page  ?orthopoly[P])

2) LegendreP(n,0) = sqrt(Pi)/(GAMMA(1/2-n/2)*GAMMA(1+n/2))
(which is interpreted as 0 if n is an odd integer).

I suspect the Maple code was entered using <maple> tags (from the maple-leaf icon).  This is now completely broken, so that it doesn't pretty-print anything, but it also doesn't let you copy anything.

Branch cut problems again.  Maple's answer using int is F(Pi/2) - F(-Pi/2) where F is the antiderivative.

> F:= unapply(int(x*sin(x)/(1+(cos(x))^4),x),x);

... a rather complicated function involving ln((r - exp(I*x))/r) where r runs over the roots
of 22 _Z^4 + _Z^8 + 4 _Z^6 + 4 _Z^2 + 1.  This will have discontinuities when
exp(I*x) = (1+t)*r, -infinity < t < 0, at approximately x = (+/-) 0.9989374565.

Export your animation as a .gif file.  The result is an animated .gif that can be viewed in any web browser.
 

Actually, you'd need to call

>  Statistics:-Sample(Statistics:-RandomVariable(...));

Otherwise you just get 0's.  I'll report this as a bug (there should at least be an error message).

Then:

> Statistics:-Sample(Statistics:-RandomVariable(Poisson(10^7)),1);

Error, (in Statistics:-Sample) NE_REAL_ARG_GT:
  On entry, rlamda must not be greater than 1.0e6: rlamda = 1e+007.

MmaTranslator[Mma][BitAnd](a,b)

The error message from pdsolve does make sense.  It means that you are not using the correct syntax for that command.  The PDE system and initial/boundary conditions must be given as sets, not just expression sequences.  Also, if you want numerical solutions you must tell that to pdsolve (I hope you're not expecting closed-form solutions for this: in any case, pdsolve without the numeric option can't handle initial and boundary conditions).  Look at the examples on the ?pdsolve,numeric help page.  I note that you defined a function Cp_air, but used Cp_Air in the definition of Prandtl.  Also, your BC4 uses a parameter T_b that has not been defined: should that be T_Bulk?

With those changes, I tried

> ans:= pdsolve({PDEsys}, {IBCsys},numeric,time=t,range=0.. R_p):

Unfortunately, I still get an error message:

Error, (in pdsolve/numeric/match_PDEs_BCs) unable to associate boundary conditions and dependent variables, system is too complex

It seems that pdsolve doesn't like your boundary conditions at x = R_p, where both D[1](T) and D[1](X) depend on both T and X.   In fact, it wouldn't even accept  D[1](X)(R_p,t)=T(R_p,t) and D[1](T)(R_p,t)=X(R_p,t), although it would accept D[1](X)(R_p,t) = X(R_p,t) and D[1](T)(R_p,t)=T(R_p,t).  I don't know if there's a work-around for this.

C'est difficile si proche au point de discontinuité, mais:

> Digits:= 20:
   Sol := dsolve({ics, sys}, 
        {v(t), z(t), m(t), vol(t), H(t), u(t)}, numeric, 
        output = listprocedure);   
   Hsol := eval(H(t), Sol);  
   Digits:= 14;
   bisect:= proc(a,b)
    local c,v;
    c:= (a+b)/2;    
    v := traperror(Hsol(c));
    if v = lasterror then bisect(a,c)
     elif v > 0.001 then bisect(c,b)
     elif c - a < 0.0001 then return a .. c
     else bisect(a,c)
    end if
   end proc; 
   bisect(0, 0.1, 10);
   fsolve(Hsol(t) = 0.001, t = %);

0.038866822386169

`#msub(mi("V"),mi("BE"))`;

Caution: There is no way this is going to be numerically stable: there is no continuous square root function on invertible matrices.  Moreover, if 0 is an eigenvalue with algebraic multiplicity >= 2, the matrix might not have a square root.  Generically there will be 8 square roots, but e.g. the identity matrix has infinitely many.

"We know it's real": yes, the matrix may be real, but the eigenvalues might not be real.  And the square root might not be real either, e.g. if the determinant of the matrix is negative there certainly can't be a real square root. 

With those cautions, here is the outline of my algorithm.

1) Find the eigenvalues r1, r2, r3 of M.  Don't use Maple's symbolic methods for this, use numerical methods.

2) Take square roots s1, s2, s3 of r1, r2, r3 respectively.  Try to make sure s1+s2, s1+s3 and s2+s3 are kept away from 0, to avoid trouble in the next step.

3) A square root of M is then
-1/((s1+s3)*(s2+s3)*(s1+s2)) * M^2
   + (s1^2+s2^2+s3^2+s1*s2+s1*s3+s2*s3)/((s1+s2)*(s1+s3)*(s2+s3)) * M
   + s1*s2*s3*(s1+s2+s3)/((s1+s2)*(s1+s3)*(s2+s3))

Use textplot to put labels at (or near) some points, or pointplot to plot some symbols there.  For example:

>  X, Y := cos(3*Pi*t), sin(2*Pi*t);
   with(plots):
   curve := plot([X, Y, t = 0 .. 1], colour=black):
   ends := pointplot(eval([X,Y], t = 0), symbol=circle, 
               colour=green),
           pointplot(eval([X,Y], t = 1), symbol = circle, 
                colour=red):
   display([curve, ends]);

Or if you want to get really fancy, insert some embedded components: a Plot (say Plot0) and a Slider (Slider0).   Set the Slider's "component properties": check "Update Continuously while Dragging", and edit "Action When Value Changes" to:

use DocumentTools, plots in
   Do ( %Plot0 = 
        display([ plot([X, Y, t = 0 .. 1], colour = black),
        pointplot(eval([X, Y], t = %Slider0/100), 
           colour = red, symbol=circle, symbolsize = 20)], 
        title = sprintf("t = %1.2f", %Slider0/100)))
 end use;

Then watch how the point moves as you move the slider.
 

M^(1/2) for a "generic" 3 x 3 matrix is going to be a mess.  But you can do it this way.
Let P(t) be the characteristic polynomial of M (its coefficients are polynomials in the entries of M).  Its roots r1, r2, r3 are the eigenvalues of M, which I'll assume are distinct.  Then M^(1/2) = Q(M) where Q is a quadratic polynomial that interpolates

Q(r1) = r1^(1/2), Q(r2) = r2^(1/2), Q(r3) = r3^(1/2)

Namely, Q(z) = r1^(1/2)*(z - r2)*(z-r3)/(r1-r2)/(r1-r3)  + r2^(1/2)*(z-r1)*(z-r3)/(r2-r1)/(r2-r3)
+ r3^(1/2)*(z-r1)*(z-r2)/(r3-r1)/(r3-r2)

If rj^(1/2) = sj (of course you have to decide which branch of the square root to use),
I get

Q(z) = -1/((s1+s3)*(s2+s3)*(s1+s2)) * z^2
   + (s1^2+s2^2+s3^2+s1*s2+s1*s3+s2*s3)/((s1+s2)*(s1+s3)*(s2+s3)) * z
   + s1*s2*s3*(s1+s2+s3)/((s1+s2)*(s1+s3)*(s2+s3))
Now all you need is to be able to solve a cubic polynomial.

Thanks, Doug.  This is very helpful.

Yes, it's amazing what a mess can happen when someone uses exact rationals instead of floats, especially with calculations in a loop.

Rather than just trying to copy your homework (or is it a take-home test?), why don't you tell us what you've done on a specific question, and where you are stuck?  We'll probably have some suggestions or hints, but (I hope) nobody here will do your work for you.