Replace dsolve by pdsolve and you get 3 solutions:

sol3[1] and sol3[2] don't depend on z.

limit(rhs(sol3[1]),rho=0,right); #Result 0
limit(rhs(sol3[1]),rho=infinity); # Result 1
## So sol3[1] satisfies all requirements regardless of the values of _C1 and _C2. WRONG:

CORRECTION: You wanted the latter limit to be zero and the first to be finite. So sol3[1] is no good.

limit(rhs(sol3[2]),rho=0,right); # Result infinity, so doesn't satisfy all requirements.

sol3[3] is considerably more complicated. It is dependent on z and involves a RootOf expression.
indets(rhs(sol3[3]),name);
indets(rhs(sol3[3]),RootOf);

You have a few problems:

1. Since laplace(T(z, t), t, s) is NOT a function of (z,t), but of (z,s), the equation sol2 should be
sol2 := subs([laplace(T(z, t), t, s) = U(z, s), T(z, 0) = sin(J*Pi*z), (D[2](T))(z, 0) = 0], sol);

2. Your boundary conditions T(0, t) = 0, T(1, t) = 0 translate into U(0,s)=0, U(1,s)=0, not what you got.

3. In order for dsolve to solve sol2 with those boundary conditions first do
   ode:=subs(U(z,s)=u(z),sol2);
   # and then
   sol3:=dsolve({ode,u(0)=0,u(1)=0});
           

4. Then you can try inverting by invlaplace:

invlaplace(rhs(sol3),s,t);  #No success.

## Comment. If you meant the boundary conditions to be D[1](T)(0,t)=0 and D[1](T)(1,t)=0, then you should do
sol3d:=dsolve({ode,D(u)(0)=0,D(u)(1)=0});
That is not, however, easier to invert.

In A the quantity l*b*c*(j+k) is an operand (or sub-operand); it is inside the square root. So here subs works.
In B the quantity j*alpha^(1/2) is not an operand (or sub-operand) of the expression 2*j*sqrt(alpha)/((j+k)*l). Thus subs won't work as intended.
Since alpha^(1/2) is an operand you can do:
C:=subs(alpha^(1/2) = beta/j, B);

You may try this, which is not terribly accurate, but seems OK:

restart;
Digits:=15:
tD := 20;
u:=4/Pi^2*( 1 - exp(-x^2*tD) ) / ( x^3*( BesselJ(0,x)^2 + BesselY(0,x)^2 ) );
series(u,x,3);
u1:=convert(%,polynom);
evalf(Int(u-u1,x=0..1/2,epsilon=1e-6));
evalf(Int(u1,x=0..1/2,epsilon=1e-5));
res1:=%%+%;
res2:=evalf(Int(u,x=1/2..infinity));
res1+res2; # Result 12.2865071846401

The integral int(u1,x=0..1/2) is finite since int(1/x/ln(x)^2,x=0..1/2) is finite.

If you change I to i (or something else) in the definition of the function pD everything works.
I get
           pDimless := 1.324262899
I is the imaginary unit.
The error message will occur in this session:

restart;

add( 5*I, I=0..3);
     
Error, illegal use of an object as a name

The variable in add (here we try I) has to be a name. I is not a name:
From the help page for I:
"I is implemented as Complex(1), and therefore, unlike many other Maple constants, type(I, name) returns false."

###### Alternatives:
1. Start the worksheet with
   restart;
   local I;
and continue without other canges.

2. Start the worksheet with
   restart;
   interface(imaginaryunit=j); #or some other name than I.
and continue without other canges.

If you do want to do the textbook steps and not just use int (lower case i) then you can do:

A:=Int(x*sqrt(2*x^4+3),x);
IntegrationTools:-Change(A,u = sqrt(2)*x^2);
value(%);
eval(%,u = sqrt(2)*x^2);
expand(%);
combine(%);
## Check with int:
value(A); #Replaces Int with int

Just a comment on your use of the (to me) confusing use of lower limit in your loop.
You have
i:=1: N:=5000:
for i from i to N
etc.

i.e. the lower limit is not written as one but i.
Both of the limits in a loop as well as the step, however, are set when the loop starts.
Take a look at this example:
restart;
j:=1:
step:=1/2:
N:=5:
for i from j to N by step do
  j:=3;
  step:=1/3;
  N:=2;
  print(i)
end do:
step,N; #They did indeed change, but that has no effect in the loop "head".
### Notice that the i-values printed are exactly the same as those in this loop:
restart;
for i from 1 to 5 by 1/2 do
  j:=3;
  step:=1/3;
  print(i)
end do:
#####
#### Now changing j to i in two examples:
## Example 1. Here it makes no difference:
restart;
i:=1:
step:=1/2:
N:=5:
for i from i to N by step do
  j:=3; #still j
  step:=1/3;
  N:=2;
  print(i)
end do:
############
## Example 2. Here we change the loop variable inside the loop to 6:
restart;
i:=1:
step:=1/2:
N:=5:
for i from i to N by step do
  i:=6; # If we still had  i:=3 the loop variable i would never get to N=5: An infinite loop.
  step:=1/3;
  N:=2;
  print(i)
end do:

LinearSolve is a member of the package LinearAlgebra.

So either use the long form LinearAlgebra:-LinearSolve or start with:

with(LinearAlgebra); # Use colon if you don't want to see the content of LinearAlgebra

pds := pdsolve(PDE, IBC, numeric, timestep = 0.01,spacestep=.5);
pds:-plot(t=0..5,x=.2);
pds:-plot(t=0..5,x=5);
#As you see here, for x>2 the solution is basically the same function of t:
pds:-animate(t=0..5);

You have 3 odes and just one unknown function g(t). Why do you expect a solution?

Secondly, if you try solving one at a time, you will find that only the last one can be solved by dsolve:

restart;  
rho:=10: #Make it concrete just for simplicity:
sys:={8*g(t)^3*diff(g(t),t$2)+4*(g(t)*diff(g(t),t))^2+1=0,rho=-1/g(t)-2*(diff(g(t),t)+t*diff(g(t),t$2))-t/(2*g(t)^3),rho=(-t/g(t))*(diff(g(t),t))^2+t/(4*g(t)^3)}
dsolve(sys[1]); #Output NULL
dsolve(sys[2]); #Output NULL
dsolve(sys[3]);
nops([%]); #6 solutions

Your procedure uses the deprecated package linalg. But that is not a problem in itself.
The commands transpose and multiply need to be replaced by
linalg[transpose] and linalg[multiply], respectively.

Alternatively, outside the procedure you can start with loading the package:
with(linalg);

I'm pretty sure that you had to do that in Maple V, Release 5, too.

The first mentioned method is the recommended one, i.e. using the long form:
linalg[transpose] and linalg[multiply].

You don't give an example of a failed run, i.e. an example of a call to Cocycle that fails.

The year in the copyright will say something, I suppose:

op(3,eval(applyrule));
    `Copyright (c) Maplesoft, a division of Waterloo Maple Inc. 2005`

On the other hand I'm sure simplify existed before this copyright year:
op(3,eval(simplify));
     system, `Copyright (c) 2000 Waterloo Maple Inc. All rights reserved.`

Take a look at
u:=unapply(U(0.15036/6.022e23,0.035,200,T),T);
# and examine
evalf[15](u~(X));
# You will notice small imaginary parts and smaller if 15 is replaced by e.g. 20.
Thus we conclude that those are due to roundoff.
A simple way to get rid of them is to take the real part as in
(Re@u)~(X);
## Thus you can do:
Statistics:-NonlinearFit(Re(U(0.15036/6.022e23,d,Theta,T)),X,Y,T,initialvalues=[Theta=200, d=0.035]);

The design of assume (and Physics:-Assume) is not made to work as eval when assumptions are equalities as in your example.
It is rather interesting though that 0 does emerge as the output from the first of these:

restart;
assume(a+b=0);
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do simplify(f(a+b)) end do;
####
restart;
Physics:-Assume(a+b=0);
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do simplify(f(a+b)) end do;
#### and also from the assuming version:
restart;
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do simplify(f(a+b)) assuming a+b=0 end do;
#######
If you remove simplify in the 3 examples above the output doesn't change.
################# Use simplify/side relations instead:
restart;
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do simplify(f(a+b),{a+b=0}) end do;
#Output 0 0 1,error as expected from ln(0).
## You can catch the errors:
restart;
for f in [x->sqrt(x^2),sin,exp,ln,x->1/x] do
  try
    res:=simplify(f(a+b),{a+b=0})
  catch:
    res:=undefined;
  end try;
  res
end do;

The description in the help page:
?dsolve,numeric,bvp

makes it clear that the methods used are finite difference schemes: Notice e.g. the options initmesh and maxmesh.
Starting with a mesh over the interval and an initial approximation to the discretized problem. That initial approximation can be given in the optional argument approxsoln.

The method used is not rkf45 as that is purely for initial value problems.