In a recent version of Maple (I tried Maple 12, 18, 2016 as a sample) try:
'rand()/rand()'; # You get one
rand()/rand(); #Highly unlikely that you get one
If you try the same two commands in Maple 8 you get 1 in both cases.
## You could try the difference also:
##In Maple 8 both of these give 0:
'rand()-rand()';
rand()-rand();
## whereas in Maple 12, 18, and 2016 you only get 0 for the first.
## If you look at the code for rand in recent versions (including 12) you will see that
RandomTools:-MersenneTwister:-GenerateInteger is used.
If you look at the code for rand in Maple 8, you will find it totally different. (use showstat(rand)).
###
While the rand procedures are surely different the reason for the different behavior in Maple 8 and Maple 2016 (say) is not that the procedures are written differently.
I tried copying the procedure from Maple 8 into Maple 2016 and gave it the name rand8.
Result:
rand8()/rand8(); #doesn't return 1 in Maple 2016.
###############
Now try this:

_n:=1;
p:=proc() global _n; _n:=_n+1; _n end proc;
p()/p(); #returns 2/3 in Maple 2016, but 1 in Maple 8
'p()/p()'; #returns 1 in both versions

Thus the change from Maple 8 to more recent versions has to do with the handling of automatic simplification.
## It should be added that p()+p() returns the same as 2*p() would, i.e. a sum like a+a is automatically simplified to 2*a.

Numerical solution has been available for quite a while.
Example with L=2:
 

restart;
ode:=diff(y(x),x$2)+lambda*y(x)=0;
bc:=y(0)=0,y(L)=0;
# Taking as an example L=2:
res:=dsolve(eval({ode,bc,D(y)(0)=1},L=2),numeric);
res(0); # You see lambda
plots:-odeplot(res,[x,y(x)],0..2);

The result of the command Sol2 := pdsolve({IBC2s, PDE2});  is NULL, i.e. Maple can't do it.
I'm using Maple 2016.2.
Try
evalb(Sol2=NULL); # true
So you are feeding build nothing, as in
build();
###
Try without IBC2s:

Sol3 := pdsolve(PDE2);
build(Sol3);

I tried your code in Maple 2016.2, Maple 2015.2, Maple 18, and Maple15:
 

restart;
with(Statistics):
X := RandomVariable(StudentT(nu));
PDF(X,0.5);

and got:

.5641895835*GAMMA(.5000000000*nu+.5000000000)/(sqrt(nu)*GAMMA(.5000000000*nu)*(1.+.25/nu)^(.5000000000*nu+.5000000000))

So I have no idea what is going on.

You will surely get some better answer, but here is a very crude way:
Insert a computation that takes a very long time at the place desired.
Interrupt the computation by using the button provided for that.
Example:
 

restart;
a:=8;
b:=89;
##
for i to 10^16 do evalf(sin(i)) end do:
##
c:=56;
d:=324;

Well, I got an answer, not great, but the supposedly good values given in your question don't seem to better at all.
 

restart;
ode_sub := diff(S(t), t) = -k1*S(t)-S(t)/T1_s;
ode_P1 := diff(P1(t), t) = 2*k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;
#ode_P1 := diff(P1(t), t) = k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;
ode_P2 := diff(P2(t), t) = -k2*(-keq*P1(t)+P2(t))/keq-k4*P2(t)-P2(t)/T1_p2;
ode_P2e := diff(P2_e(t), t) = k4*P2(t)-P2_e(t)/T1_p2_e;
ode_system := ode_sub, ode_P1, ode_P2, ode_P2e;
s0 := 96304.74567; k2 := 10^5; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1;
init := S(0) = s0, P1(0) = 0, P2(0) = 0, P2_e(0) = 0;
ode_system;
T := [0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58];
####
data_s := [96304.74567, 77385.03700, 62621.83067, 51239.94333, 42663.82367, 35084.74100, 28480.28367, 23066.01467, 18774.73700, 15179.13700, 12278.50767, 9937.652000, 8046.848333, 6521.242000, 5287.811667, 4277.779000, 3466.518333, 2835.467000, 2297.796333, 1861.249667, 1529.654000, 1235.353000, 999.6626667, 826.2343333, 667.9480000, 559.9230000, 449.2790000, 376.4860000, 289.1203333, 236.1483333];
####
data_p1 := [0.86e-1, 3.904, 26.975, 31.719, 41.067, 46.779, 52.115, 43.101, 44.344, 41.094, 36.523, 27.742, 26.543, 28.062, 22.178, 21.303, 14.951, 17.871, 11.422, 12.051, 9.232, 6.817, 6.1, .717, 1.215, 6.146, .772, .375, 2.595, .518];
####
data_p2_p2e := [-3.024, 22.238, 61.731, 103.816, 132.695, 159.069, 167.302, 160.188, 158.398, 152.943, 146.745, 135.22, 132.145, 120.413, 107.864, 95.339, 90.775, 81.828, 71.065, 70.475, 62.872, 49.955, 40.858, 42.938, 41.311, 35.583, 31.573, 29.841, 29.558, 21.762];
#### The data are vastly different in size, so it might be good to put weights on the residuals.
#### I have done that in the uncommented version of resid in QM
max(data_s);
max(data_p1);
max(data_p2_p2e);
####
beta:=10: Tr:=2:
K:=unapply(evalf(cos((1/180)*beta*Pi))^(t/Tr),t);
KT:=<K~(T)>;
####
resM:=(T1_p1xx, k1xx, keqxx, k4xx)->
  dsolve(eval({ode_system,init},{T1_p1=T1_p1xx, k1=k1xx, keq=keqxx, k4=k4xx}),numeric,_rest);
#### Dividing the data elementwise by the elements in KT:
P1d:=<data_p1>/~KT;
P2_P2ed:=<data_p2_p2e>/~KT;
Sd:=<data_s>/~KT;
###############################################################################
QM:=proc(T1_p1, k1, keq, k4) option remember; local res,A,P1v,P2v,P2e_v,Sv,resid;
    userinfo(1,procname,printf("Trying %a\n",[_passed]));
    try
     res:=resM(T1_p1, k1, keq, k4,maxfun=10^8,output=Array(T),stiff=true); #stiff used
     A:=res[2,1];
     ## 
     if A[-1,1]=`?` then error "cannot evaluate the solution" end if;
     catch "cannot evaluate the solution":
       return [10^6$(3*nops(T))];
     end try;
     P1v:=A[..,2];
     P2v:=A[..,3];
     P2e_v:=A[..,4];
     Sv:=A[..,5];
     #resid:=[P1v-P1d,P2v+P2e_v-P2_P2ed,Sv-Sd];# No weights
     resid:=[(P1v-P1d)/max(P1d),(P2v+P2e_v-P2_P2ed)/max(P2_P2ed),(Sv-Sd)/max(Sd)];#Weight put on all
     [seq(seq(resid[i][j],i=1..3),j=1..nops(T))]
     end proc:
q:=[seq(subs(_nn=n,(proc(T1_p1, k1, keq, k4) QM(_passed)[_nn] end proc)),n=1..3*nops(T))]: 
#################################################################################
## The wrapping of Codetools:-Usage around the LSSolve command is done for measuring time etc.
##
forget(QM):
sol:=CodeTools:-Usage(Optimization:-LSSolve(q,initialpoint=[0.03, 0.02, .5, .5],output=solutionmodule));
sol:-Results();
sol:-Settings();
Number of calls to QM:
nops([indices(op(4,eval(QM)))]);
params:=convert(sol:-Results(solutionpoint),list);
RES:=resM(op(params),stiff=true);
pS:=plot(T,data_s,style=point,symbolsize=20,symbol=solidcircle):
pP1:=plot(T,data_p1,style=point,symbolsize=20,symbol=solidcircle):
pP22e:=plot(T,data_p2_p2e,style=point,symbolsize=20,symbol=solidcircle):
plots:-display(pS,plots:-odeplot(RES,[t,S(t)],0..T[-1]),labels=[t,S]);
plots:-display(pP1,plots:-odeplot(RES,[t,P1(t)],0..T[-1]),labels=[t,P1],thickness=3); 
plots:-display(pP22e,plots:-odeplot(RES,[t,P2(t)+P2_e(t)],0..T[-1]),labels=[t,P2+P2_e]);
The values given as "the solution" don't fit very well:
RESsol:=resM(36.8,0.000438,2.7385,0.0385,stiff=true);
plots:-display(pS,plots:-odeplot(RESsol,[t,S(t)],0..T[-1]),labels=[t,S]);
plots:-display(pP1,plots:-odeplot(RESsol,[t,P1(t)],0..T[-1]),labels=[t,P1],thickness=3); 
plots:-display(pP22e,plots:-odeplot(RESsol,[t,P2(t)+P2_e(t)],0..T[-1]),labels=[t,P2+P2_e])

As the setting is above the result is
params := [0.0501239039106269, 0.0115045966178040, 0.663492002313513, 0.434185918627229]
Remember that the order is [T1_p1, k1, keq, k4].
The graphs are not great. My only comment to that is one I have given before:
If the data are real data, then it is unlikely that your model is all that good. After all it is linear, so without knowing anything about the origin of the model, my guess is that it is used because it is the simplest one could come up with. Never a bad start though.

If weights are not put on the residuals the graph for S looks almost perfect, but that is no miracle since data are so big.

First of all you can only work on a finite interval. Thus you must replace infinity by some (in the context) large number.
I have used subs(infinity=4,{bc}) below, but you can just remove infinity and write 4 instead.
Secondly, the second derivative in bc should be written (D@D)(f)(0) or alternatively as (D@@2)(f)(0). The latter version is easier to use for higher derivatives.
Thirdly, cosine is spelled cos, not Cos.
Fourthly, use a midpoint method.
After these corrections the following works. Notice that I use infinity=4.

res:=dsolve({eqn1,eqn2,eqn3} union subs(infinity=4,{bc}),numeric,method=bvp[midrich]);
plots:-odeplot(res,[eta,f(eta)]);

Try this:
 

f:=x->piecewise(0<x and x<0.5,x,0.5<x and x<1,1-x, 0);
g:=unapply(convert(f(x),Heaviside),x);
# Then compare:
diff(f(x),x);
diff(g(x),x);

Just use evalc as in
v := evalc(Vector(4, [Re(V[1, 1]), Im(V[1, 2]), Re(V[2, 1]), Im(V[1, 1])]));

evalc assumes itself that variables are real.
No need for making assumptions yourself.
You had a syntax error: When I copy the assume line I get
`assume `*{psi(t), 'real'}
which is utter nonsense:
(1) Curly braces are used for defining sets only.
(2) The appearance of the multiplication sign is due to the silly convention in 2D math input (the unfortunate default) that a space means multiplication.
But the good news: You don't need assume here.

You could normalize the rows like this:

M := Matrix([[1,2],[3,4]]);
<seq(M[i]/LinearAlgebra:-Norm(M[i],2),i=1..2)>;

Ij you have an initial value problem you can set abserr and relerr.
?dsolve,numeric,ivp  and the reference in that page to Error Control for Numerical Solution of IVPs in Maple
If you have a boundary value problem you can set abserr.
?dsolve,numeric,bvp

If you are are using (or in effect using) conjugate for getting the bars over symbols your problem is due to the fact that gamma is Euler's constant in Maple, thus it is real. Try
evalf(gamma);
You can get over that by starting your worksheet with a local declaration like this:
restart;
local gamma;
##And then proceed as in
conjugate(gamma);

Try in the symbolic case to replace V^%T by V^(-1), which is the correct matrix in the general case.
Then all should be just fine, i.e.
V . DiagonalMatrix(E) . V^(-1);
returns Omega exactly.
###################
A comment.
The procedure IsMatrixShape is somewhat misleading.
Consider this simple example where M in fact happens to be symmetric, but doesn't have that shape as one of its options. Recall that eigenvalues of real symmetric matrices are real.
 

M:=Matrix(3,(i,j)->i+j);
IsMatrixShape(M,symmetric);
MatrixOptions(M); # shape=[ ]
Ms:=Matrix(M,shape=symmetric);
MatrixOptions(Ms); # shape = [symmetric]
LinearAlgebra:-Eigenvalues(evalf(M));
LinearAlgebra:-Eigenvalues(evalf(Ms));

I think that the simplest way to achieve that is by going straight to the definition of the taylor expansion:
 

f1 := taylor(tanh(sqrt(q)*z/y+x*m/y), m, 2);
f2:=convert(f1,polynom);
f:=m->tanh(sqrt(q)*z/y+x*m/y);
f3:=add((D@@i)(f)(0)*m^i/i!,i=0..1);#The definition
simplify(convert(f2-f3,exp)); # 0

What you have is a text string. Nothing whatsoever inside the string is evaluated.
So even if n:=2, the string "n" remains the string "n".
You could use typest as in

plot(A[1](t), t=0..tf, labels = ["Time t", typeset( n, " times derivative of ", X[1](t) ) ] );

And since I notice that you don't want X[1](t) evaluated, you could do:

plot(A[1](t),t=0..tf,labels=["Time t",typeset(n," times derivative of ", 'X[1](t)')]);

Notice the unevaluation quotes '   '.
#Comment: By 2 times derivative you really mean the second derivative, I'm sure.