You can use normal, expanded instead, either directly or via applyop:

normal(m,expanded);
applyop(normal,3,m,expanded);

There are some syntactical problems, among them missing multiplication signs.

You have 6 boundary conditions instead of the expected 5 (the sum of the orders 3+2). Does that mean that lambda has to be determined?
Furthermore you have
theta(0)=0, theta(0)=1+sigma*D(theta)(0)
thus if sigma<>0 then theta and theta' are both zero at 0, so theta(eta)=0 will be the trivial solution.

The system as you have it after corrections:
sys:=[diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta,eta)+(2*n)/(n+1)*(1-diff(f(eta),eta)^2) = 0, 1/(Pr) *diff(theta(eta),eta,eta)+f(eta)*diff(theta(eta),eta)-(2*p)/(n+1)*diff(f(eta),eta)*theta(eta) = 0, f(0) = 0, D(f)(0) = 1+lambda*(D@@2)(f)(0), D(f)(10) = 1, theta(0)=0, theta(0)=1+sigma*D(theta)(0),theta(10)=1];

###

When successful and if you want f''(eta) do:
plots[odeplot](sol1, [eta, diff(f(eta),eta,eta)], color = red,axes=boxed);

For the present case this would do it:

remove(has,eq2,cos);

# The following is more general:

remove(depends,eq2,x);

Your code reminds me strongly of a webinar by Sam Dao available from Maplesoft: Recorded webinars.

You have the problem that what you define as ST contains an x that shouldn't appear. The space variable is supposed to be discretized.

You have:

ST := k -> (u[k+1](t)-2*u[k](t)+u[k-1](t))/h^2+2*u[k](t)-3*k*h-2*(diff(u[k](x, t), t));

The last term should be reconsidered.

Another problem: You assign to t[max] in the beginning. That makes t a table implicitly. This causes a problem at the end.
In general one should avoid using a name (like t) and its indexed version (as in t[max] ) in the same worksheet.
Change all occurrences of t[max] to something else, e.g. tmax or (t__max).

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.