You should begin by replacing all square brackets i.e. [ ], by parentheses, ( ), since square brackets are used for other things e.g. for making lists:  [a,b,c] is a list and so is [a] .

Secondly I see missing multiplication signs (*) in e.g. tau_mafw, tau_mafo, tau_mafg.

You cannot use a name like S_wfi both as a function name as it is in S_wfi(x,t) and also as a constant as it is in e.g. IBC.
Maybe the appearance of S_wfi(x,t) is an error and should have been S_wf(x,t) ?

After those corrections (I included the replacement of S_wfi(x,t) with S_wf(x,t) ) the variables can be found by these 3 steps

indets({Eq1, Eq10, Eq11, Eq12, Eq2, Eq3, Eq4, Eq5, Eq6, Eq7, Eq8, Eq9},specfunc(diff));
op~(1,%);
vars:=% minus indets(%,specfunc(diff));

But you will have no success with
pdsolve({Eq1, Eq10, Eq11, Eq12, Eq2, Eq3, Eq4, Eq5, Eq6, Eq7, Eq8, Eq9, IBC},vars);

which doesn't surprise me at all.
If you leave out IBC then pdsolve will work for a long time, but I interrupted the computation since I didn't expect an answer anyway.

Down to basics I would do:

restart;
v1:=<1,0,0>: v2:=<0,1,0>: v3:=<0,0,1>:
u1:=<0,1,0>: u2:=<0,0,1>:
A:=<v1|v2|v3>;
#We need to find out if both of the equations A.x = u1 and A.x = u2 have solutions.
# That is the same as asking if  A.X = <u1|u2> has a solution X.
# Thus we form the augmented matrix:
T:=<A|u1|u2>;
## and then simply:
LinearAlgebra:-Rank(T) - LinearAlgebra:-Rank(A); #Yes if 0, No if greater than zero
## You could solve, but that is not necessary here.
#LinearAlgebra:-LinearSolve(A,<u1|u2>);

I have no problem in Windows 10 and with the 64 bit version of Maple 2016.1.

In Maple 2016 the two prints in each procedure produce the same, thus you get
                              0, 0
                              1, 1
                              0, 0
                              1, 1
In Maple 8 it is as you report for Maple V Release 4:

                                 0, 0
                                 0, 1
                                 0, 0
                                 1, 1
So it seems that assign treats environtment variables different than other variables in early versions of Maple.
The help page for Environment Variables in Maple V Release 4 (and in Maple 8) starts by saying:
     "Environment variables can be used in a simple assignment."

But the help page for Maple 2016 says the exact same thing.
In Maple V Release 4 and in Maple 8 the global version
_EnvX:=0; assign('_EnvX=1');
doesn't change _EnvX to 1 either as does Maple 2016.

In Maple 2016 all four produce the same output:

map((a::uneval,b)->'args',[a,b,c,d],1..4,x);
map((a,b::uneval)->'args',[a,b,c,d],1..4,x);           
map((a,b)->'args',[a,b,c,d],1..4,x);
map((a,b)->args,[a,b,c,d],1..4,x);

## The output is
      [a, 1 .. 4, x, b, 1 .. 4, x, c, 1 .. 4, x, d, 1 .. 4, x]

When you say the product is Maple V, which release is it?

The output is the same in Maple 12 as in Maple 2016, but Maple 8 returns ['a, 1 .. 4, x', 'b, 1 .. 4, x', 'c, 1 .. 4, x', 'd, 1 .. 4, x'] in the first case, and [a, 1 .. 4, x, b, 1 .. 4, x, c, 1 .. 4, x, d, 1 .. 4, x] in the other three.

The function ux has complex values, so plot can only plot the real and/or the imaginary parts:

plot([Re@ux,Im@ux], 0..10, labels=[x,"Re(u),Im(u)"],color=[green,red],caption="t=0.8");

I suppose that the Digits setting at 5 is a misprint and ought to have been at least 50 as your data appears to have that many digits.
At Digits:=50 I don't get anything like the graph you show, not for Re(u) nor for Im(u).
Raising Digits to 60 still doesn't produce that graph. There are huge values at x > about 7.
Raising Digits to 100 produced this:



The numerical inversion of the laplace transform is notoriously difficult, I understand.

nb/lhs(hh)*rhs(hh);

I never use patmatch, so I cannot be of much help there. Notice though that you don't have a square root at all.

indets(f1,radical); #A fourth root
nops(f1); # 3
whattype(f1); # `*` a product
op(f1); #The 3 operands in the product
op(2,f1); # We can extract the radical and the polynomial part from this factor:
rad,pol:=selectremove(type,op(2,f1),radical);
rad;
pol;
## You shouldn't forget the two other factors, though:
op(1,f1);
op(3,f1);

There seems to be a solve problem as well as a LinearSolve problem. The reason is not clear to me right now.
(See note 3 at the bottom for another simple way out).
But your system SYS2 with the boundary conditions bcs2 and bcs22 can be solved like this.

Define SYS2, bcs2, and bcs22 as you did.
Then do:

Hdiff:={diff(h1(theta),theta$4),diff(h2(theta),theta$3),diff(h3(theta),theta$4)};
#solve(SYS2,Hdiff); #Does it finish?
A,b:=LinearAlgebra:-GenerateMatrix(convert(SYS2,list),[op(Hdiff)]);
# Try
LinearAlgebra:-Determinant(A); # Pretty small
LinearAlgebra:-ConditionNumber(A); #Pretty large
#LinearAlgebra:-LinearSolve(A,b); #Does it finish?
## That could be the reason for solve and fsolve not to return in a reasonable time; see note at bottom, though.
## Instead of using LinearSolve, we just do
sol:=A^(-1).b;
sys:=[op(Hdiff)]=~convert(sol,list):
res:=dsolve({op(sys)} union bcs2 union {bcs22},numeric);
plots:-display(Array([seq(plots:-odeplot(res,[theta,h(theta)]),h=[h1,h2,h3])]));
###############
Note added. The problem for LinearSolve seems to have something to do with the vector b rather than with the matrix A:
This causes absolutely no problem:
LinearAlgebra:-LinearSolve(A,<x,y,z>);
sol:=eval(%,{x=b[1],y=b[2],z=b[3]});
#Then continue as before.
Note 2 added: I tried some freezing when using solve with not that much success:
Hsbs:=Hdiff=~freeze~(Hdiff);
Fus:=indets(SYS2,function) minus Hdiff;
Fsbs:=Fus=~freeze~(Fus);
SYS2F:=subs(Hsbs,Fsbs,SYS2);
sol:=solve(SYS2F,rhs~(Hsbs)): #Use colon since it is huge.
sys:=thaw(sol);
dsolve(sys union bcs2 union {bcs22},numeric); #Complains about singularity at left endpoint
dsolve(sys union bcs2 union {bcs22},numeric,method=bvp[midrich]); #Complains about singularity
## I wouldn't trust these last two "results".
###############################
Note 3: Freezing expressions of type `^`  helps greatly:
SYS2FF:=evalindets(SYS2,`^`,freeze);
sys:=thaw(solve(SYS2FF,Hdiff));# FAST
# Now proceed with
res:=dsolve(sys union bcs2 union {bcs22},numeric);

You have 3 independent variables. You have appearing X(T,Z) and Y(T,R), so 3 independent variables T,Z,R.
That is one two many for pdsolve/numeric.
Did you mean Z when you wrote R?
The syntax would have been
pdsolve({EQ[1],EQ[2]}, BC[1] union BC[2],numeric);

What is Xstar?

1. cosine and sine only appear as cos(g) and sin(g) and g:=Pi/6, so its irrelevant that they appear in input.
2. I don't understand what this is doing in your question: "Failed to load the worksheet /maplenet/convert/Sc.mw . "
3. I looked at your worksheet. Since immediate success is not achieved even after using a midpoint method, as you are told to do by the error message (add method=bvp[midrich] ), and since you are using continuation, I suggest beginning more modestly.
Eq1 has only got the dependent variable f. Thus I suggest studying that first with the boundary conditions
bcsf:=f(0) = 0, D(f)(0) = lambda, D(f)(5) = 1;
#
Continuation in the parameter d will only work if dsolve can solve the problem when d=0, because it works from there and up to d=1 in steps using previously found results as approximate solutions along the way. The problem you really want to solve corresponds to d=1.
So begin by trying d = 0 and one of your lambda values (don't assign to lambda, use eval).
# You can try
res:=dsolve(eval({Eq1,bcsf},{lambda=-2.6,d=0}),numeric,method=bvp[midrich]);
# but you get an error saying

Error, (in dsolve/numeric/bvp) unable to achieve requested accuracy of 0.1e-5 with maximum 128 point mesh (was able to get 0.68e-2), consider increasing `maxmesh` or using larger `abserr`

This, however, may very well lead to a wild goose chase: increasing maxmesh without success. Increasing abserr all the way up to 0.04 works, though:
res:=dsolve(eval({Eq1,bcsf},{lambda=-2.6,d=0}),numeric,method=bvp[midrich],abserr=.04);
##
There is a long way from that to d=1, i.e. using continuation. No success here:
res:=dsolve(eval({Eq1,bcsf},{lambda=-2.6}),numeric,method=bvp[midrich],abserr=.04,continuation=d);
## One may wonder if continuation in the coefficient of eta*f'''(eta) really is a good idea.
##
I tried with (in my view) the absurdly large absolute error abserr=3 (and on the real problem, i.e. d=1):
res1:=dsolve(eval({Eq1,bcsf},{lambda=-2.6,d=1}),numeric,method=bvp[midrich],abserr=3);
#Success, but is it? What I got was
plots:-odeplot(res1,[seq([eta,diff(f(eta),[eta$i])],i=0..2)],legend=[f,seq(Diff(f,[eta$i]),i=1..2)]);


####### Pursuing this further still with abserr=3, no continuation, and lambda=-2.6:
resfull:=dsolve(eval(dsys,{lambda=-2.6,d=1}),numeric,method=bvp[midrich],abserr=3);
plots:-odeplot(resfull,[[eta,f(eta)],[eta,phi(eta)],[eta,theta(eta)]]);


This works for the two other values of lambda as well, but it does bother me that I have to use abserr=3.

Maybe frontend is what you need:

eq := (diff(g(z), x))*g(z)^3+(diff(g(z), z, z))*g(z)^4+(diff(g(z), z, z, z))*g(z)^5+(diff(g(z), z))/g(z)^2;
frontend(coeff,[eq,g(z)^4]);
frontend(coeff,[eq,g(z)^(-2)]);

Although equality assumptions are allowed (I think primarily for reasons of consistency), they don't work as evaluation.
Simple example:
restart;
sqrt(a^2) assuming a=2;

The output is a not 2.
Thus the properties the integer 2 has (like being positive) have been used, but a has not been replaced by 2.

Set Digits higher than the default whish is 10, doing like this: Digits:=50: (or whatever you like!).
Alternatively do
restart;
evalf[50](BesselI(0,31.7)); # or use another integer than 50, e.g. 15
## I get
                 4.1615582165614032637527525720570080615945783690515*10^12
# Notice 10^12, not 10^20.
# Similarly
evalf[50](BesselI(0,60));
# gives me  
       5.8940770556098011682788174403339047379789830203162*10^24

First a theoretical remark:
Rewrite your ode as a system:
DEtools[convertsys](ode,{},theta,t,Y,YP);
   [[YP[1] = Y[2], YP[2] = sin(Y[1])*cos(Y[1])-9.8100*sin(Y[1])], [Y[1] = theta(t), Y[2] = diff(theta(t), t)], undefined, []]
You have an autonomous system in Y= (Y[1],Y[2]). The system has an equilibrium at Y = (Pi,0).
If you start at (Pi,0) you will never leave, if you don't, you will never get there!
So the short answer is: No solution!
##Further comments:
The equilibrium (Pi,0) is a saddle point, thus unstable.
The jacobian (where I use g for 9.81):
J:=unapply(VectorCalculus:-Jacobian([Y[2],sin(Y[1])*cos(Y[1])-g*sin(Y[1])],[Y[1],Y[2]]),[Y[1],Y[2]]):
J(Pi,0);
LinearAlgebra:-Eigenvalues(%); #One positive and one negative, so a saddle point.
A saddle point has 4 separatrices corresponding to solutions satisfying Y(t) -> (Pi,0) as t -> infinity or as t->-infinity.
The directions of approach near (0,Pi) are given by the eigenvectors.
Below we find a solution which gets very close to a separatrix. The solution is approaching (0,Pi) and gets very close, but leaves the vicinity again. See the animation.

So try this:

restart;
ode:=sin(theta(t))*cos(theta(t))-9.8100*sin(theta(t))-(diff(theta(t), t, t)) = 0;
sol:=dsolve({ode, theta(0) = 2*Pi*(1/3), D(theta)(0) = th1}, numeric,parameters=[th1],output=listprocedure);
TH,TH1:=op(subs(sol,[theta(t),diff(theta(t),t)]));
## The way parameters are set:
sol(parameters=[1.0341*Pi]); #Example
plots:-odeplot(sol,[[t,theta(t)-Pi],[t,diff(theta(t),t)]],t=0..5) ;
##Now make a procedure with output we want to be [0,0]:
p:=proc(th1,t0) option remember; sol(parameters=[th1]); [TH(t0)-evalf(Pi),TH1(t0)] end proc;
## Because of LSSolve (or fsolve) we need these:
p1:=proc(th1,t0) p(_passed)[1] end proc;
p2:=proc(th1,t0) p(_passed)[2] end proc;
## Testing:
p1(1.0341*Pi,1.5);
p(1.0341*Pi,1.5);
# I haven't had success with fsolve, so I use LSSolve:
res:=Optimization:-LSSolve([p1,p2],initialpoint=evalf([1.0340*Pi,1.5]));
res[2][1]/evalf(Pi);
plots:-odeplot(sol,[[t,theta(t)-Pi],[t,diff(theta(t),t)]],t=0..6) ;
## Phase plane:
plots:-odeplot(sol,[[theta(t),diff(theta(t),t)]],t=0..10,frames=50,caption="Very close to a separatrix") ;

################## Finding an equation for the separatrices:
ode:=sin(theta(t))*cos(theta(t))-9.81*sin(theta(t))-(diff(theta(t), t, t)) = 0;
## ode implies that the following is zero:
lhs(ode)*diff(theta(t),t);
#So we have
eq:=int(%,t)=C; #Where C is a constant
#Finding C corresponding to a separatrix for (Pi,0):
eval[recurse](convert(eq,D),{t=0,theta(0) = Pi, D(theta)(0)=0});
#So an equation for a separatrix is
eq2:=subs(diff(theta(t),t)=theta_p,theta(t)=theta,eval(eq,C=-9.81));
plots:-implicitplot(eq2,theta=0..15,theta_p=-7..7,gridrefine=2,labels=[theta,theta__p]);