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.

The response from the LinearSolve command is a "return unevaluated". That is, nothing happens, you just get the input back although M is evaluated.
This would happen if you forgot to execute with(LinearAlgebra).
So try again:
 

restart;
with(LinearAlgebra):
M:=<<1,2,3>|<2,2,2>|<4,1,0>|<5,4,1>>;
LinearSolve(M);

You can make sure that the package LinearAlgebra has been activated by replacing the colon by a semicolon. If the package content is shown when you use semicolon, things should be in order.
An alternative to using with(LinearAlgebra) is to use the long form:
 

restart;
#with(LinearAlgebra):
M:=<<1,2,3>|<2,2,2>|<4,1,0>|<5,4,1>>;
LinearAlgebra:-LinearSolve(M);

You already got 2 excellent answers.
Nevertheless, and for the fun of it, I shall give a third way.
 

restart;
eq := (1/4)*x*(3-2*x)^2/(1-x)^3 = 18325.73901+exp(36.58420421-10902.09286/T);
ode:=diff(subs(x=x(T),eq),T);
ic:=x(300)=fsolve(eval(eq,T=300.),x);
res:=dsolve({ode,ic},numeric);
plots:-odeplot(res,[T,x(T)],300..800,size=[800,default]);
res(500); #value of x at T=500

The person at your university who provided you with the program ought to have the purchase code or should know what to do.

In Maple try
?backup
You see that you should do this:

From the File menu, select Recent Documents, then Restore Backup. All backed up files are opened.

The following version of my answer given earlier is about 10 times faster.
It doesn't use option parameters in dsolve/numeric. It can't because it uses output=Array(T1), where T1 is the given t-data (a list).
On the data we have it takes just less than 0.5 minutes on my computer, where the parameters + listprocedure version took 4.5 minutes.
 

restart;
T1 := [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, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100]:

P1d1 := [-0.8685347292816209e-3, 66.83841848683134, 111.424044735605, 139.4041294965286, 155.062049090439, 161.7293955199091, 161.9517800941181, 157.6868656783518, 150.4175354290345, 141.2545066744467, 131.0251276136712, 120.3351733051009, 109.6105958753724, 99.15971094060681, 89.1849815138892, 79.80813011118238, 71.10663373292822, 63.11585091923417, 55.83128037844, 49.24042186072761, 43.31235428504942, 38.00371754713981, 33.27679003320984, 29.07632002918529, 25.36059655209477, 22.08457240828178, 19.20460745561158, 16.67539869671402, 14.4614396803485, 12.52248311609639, 10.83757549915796, 9.36708834933514, 8.088032224392839, 6.978649948380624, 6.015427486670388, 5.18019389659066, 4.459752143602849, 3.833458379711213, 3.296146126778823, 2.834176028037793, 2.432779628368156, 2.086515038739586, 1.788690970907987, 1.533396481460015, 1.314356306978484, 1.12730777572421, .9616785898396352, .8260162377043778, .7063555982247525, .6010919595382802, .5139462672093198]:

Sd1 := [9999.99913146527, 8202.338628704818, 6727.833761289698, 5518.396930229777, 4526.371761304796, 3712.685113725112, 3045.268751278301, 2497.831339895754, 2048.805817405232, 1680.499031190071, 1378.401697618653, 1130.613667949291, 927.3657548306696, 760.6567913452087, 623.9185918092153, 511.7575147286922, 419.7595064851229, 344.3032822013375, 282.4083274180995, 231.6408689064719, 190.0007294821676, 155.8435093977208, 127.8303140245707, 104.8488496263611, 85.99905541436985, 70.53919201906764, 57.86040776255557, 47.4590859030469, 38.92828403470473, 31.92689759464799, 26.19042899044342, 21.4822256595649, 17.62001555044496, 14.45336300140778, 11.85487894728136, 9.722608317770803, 7.975833977605903, 6.539549720337241, 5.364864331917502, 4.40300772819562, 3.611018434464852, 2.960857591104723, 2.427746882766981, 1.991319572019064, 1.633765408748718, 1.341633570447127, 1.097047090111865, .9027584989710459, .7402161952815224, .6041949051082868, .4955658248846614]:

P2_P2ed1 := [-0.17370694585632418e-2, 269.58999604287425, 454.8469765120931, 575.9475396050403, 648.7737935745041, 685.7457798251972, 696.4686231070878, 688.3821550997267, 667.1907577733688, 637.2317161484453, 601.7834093079005, 563.2932079626577, 523.544273610169, 483.8462546459006, 445.1115231227346, 407.94720742960124, 372.7612377498941, 339.7926684547662, 309.1366947159386, 280.8244319793804, 254.81459910501462, 231.0152072214917, 209.3272228584818, 189.60398892316664, 171.71774749517613, 155.52918717277277, 140.89938734410705, 127.68548032394494, 115.76486486199599, 105.00703614849681, 95.32226228958285, 86.58602900058177, 78.7098815088867, 71.61070067323783, 65.20472435079975, 59.42168320280972, 54.20395463499331, 49.481585328835266, 45.21798670134083, 41.36187311862125, 37.86178597558813, 34.68684377541566, 31.80543424501873, 29.189125200053248, 26.81023517978529, 24.646456502685155, 22.664089588518717, 20.868790019282276, 19.222201203537946, 17.711606686397563, 16.33577524343594]:
##############################################################################
ode_sub := diff(S(t), t) = -k1*S(t)-S(t)/T1_s;
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;
s0:=10000;
k2 := 1000; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1;

ode_system := ode_sub, ode_P1, ode_P2, ode_P2e;
init:=S(0)=s0,P1(0)=0,P2(0)=0,P2_e(0)=0;
###############################################################################
resM:=(T1_p1xx, k1xx, keqxx, k4xx)->
  dsolve(eval({ode_system,init},{T1_p1=T1_p1xx, k1=k1xx, keq=keqxx, k4=k4xx}),numeric,output=Array(T1),maxfun=10^6);
###############################################################################
QM:=proc(T1_p1, k1, keq, k4) option remember; local res,A,P1v,P2v,P2e_v,Sv,resid;
    printf("Trying %a\n",[_passed]);
    try
     res:=resM(T1_p1, k1, keq, k4); #This won't protest if there is a problem.
     ## So we force it to do so:
     if has(res,`?`) then error "cannot evaluate the solution" end if;
     catch "cannot evaluate the solution":
       return [10^6$(3*nops(T1))];
     end try;
     A:=res[2,1];
     P1v:=A[..,2];
     P2v:=A[..,3];
     P2e_v:=A[..,4];
     Sv:=A[..,5];
     resid:=[P1v-<P1d1>,P2v+P2e_v-<P2_P2ed1>,Sv-<Sd1>];
     [seq(seq(resid[i][j],i=1..3),j=1..nops(T1))]
end proc:
q:=[seq(subs(_nn=n,(proc(T1_p1, k1, keq, k4) QM(_passed)[_nn] end proc)),n=1..3*nops(T1))]: 
#################################################################################
## The wrapping of Codetools:-Usage around the LSSolve command is done for measuring time etc.
##
sol:=CodeTools:-Usage(Optimization:-LSSolve(q,initialpoint=[1,1,1,1]));
params:=convert(sol[2],list);
res:=resM(op(params));
p1:=plot(T1,P2_P2ed1,style=point,symbolsize=20,symbol=solidcircle);
p2:=plots:-odeplot(res,[t,P2(t)+P2_e(t)],0..T1[-1]);
plots:-display(p1,p2);

Take a look at
?solve,details
Example:
solve(mul(x-k,k=1..20),x,maxsols=3);
solve(mul(x-k,k=1..20),x,maxsols=1);

The following works now.
 

restart;
ode_sub := diff(S(t), t) = -k1*S(t)-S(t)/T1_s;
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;
s0 := 10000;#edited: was 100000
k2 := 1000; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1;
ode_system := ode_sub, ode_P1, ode_P2, ode_P2e;
indets({ode_system},name); #Just a check
init:=S(0)=s0,P1(0)=0,P2(0)=0,P2_e(0)=0;
####
res:=dsolve({ode_system,init},numeric,parameters=[T1_p1, k1, keq, k4],output=listprocedure,maxfun=10^6);
####
P1fu,P2fu,P2e_fu,Sfu:=op(subs(res,[P1(t),P2(t),P2_e(t),S(t)]));
####
T:=<0,2,4,6,8>;
Sd:=<9999.99913146527,8328.870587730016,6937.009129218748,5777.745632133724,4812.209983843559>;
P1d:=<0.0,67.86790056712294,114.88787098501874,145.95438088662502,164.85650644237887>;
P2_P2ed:=<0.0,271.68492651947497,461.9130396605823,589.3710176125417,668.9967533337124>; # data from P2(t)+P2_e(t)
####
Q:=proc(T1_p1, k1, keq, k4) option remember; local P1v,P2v,P2e_v,Sv,resid;
   res(parameters=[T1_p1, k1, keq, k4]);
   P1v:=P1fu~(T);
   P2v:=P2fu~(T);
   P2e_v:=P2e_fu~(T);
   Sv:=Sfu~(T);
   resid:=[P1v-P1d,P2v+P2e_v-P2_P2ed,Sv-Sd];
   [seq(seq(resid[i][j],i=1..3),j=1..5)]
end proc:
q:=[seq(subs(_nn=n,(proc(T1_p1, k1, keq, k4) Q(_passed)[_nn] end proc)),n=1..15)]: 
Q(1,2,3,4); #Check
q[4](1,2,3,4); #Check
### q is a list of procedures each producing a residual.
### Knowing a reasonable initialpoint (guess for the parameters) could help.
#Will find one with fsolve on 4 of the 15 q procedures

L:=fsolve(q[2..5],[10,0.02,4,4]); #To get an idea for a guess
solLS:=Optimization:-LSSolve(q,initialpoint=L);
res(parameters=convert(solLS[2],list));
with(plots):
pS:=plot(T,Sd,style=point,symbolsize=20,symbol=solidcircle);
pP1:=plot(T,P1d,style=point,symbolsize=20,symbol=solidcircle);
pP22e:=plot(T,P2_P2ed,style=point,symbolsize=20,symbol=solidcircle);
display(pS,odeplot(res,[t,S(t)],0..8),labels=[t,S]);
display(pP1,odeplot(res,[t,P1(t)],0..8),labels=[t,P1]);
display(pP22e,odeplot(res,[t,P2(t)+P2_e(t)],0..8),labels=[t,P2+P2_e]);

Here is the last plot:

When I try in my Maple 2016.2 I certainly get a result.
Full code:
 

restart; 
with(PDEtools);
U := diff_table(u(x, y, z));
pde[1] := x*y*U[z]+x*U[x]+2*y*U[y] = 0;
bc[1] := eval(U[], z = 0) = x^2+y^2;
sys[1] := [pde[1], bc[1]];
sol:=pdsolve(sys[1]); #NOTICE cubic in _Z
allvalues(sol); # Three solutions

For sol I get
sol := u(x, y, z) = (x*y-3*z)*(RootOf(_Z^3*y-x^3*y+3*x^2*z)*x^2+(x*y-3*z)*y)/(y*RootOf(_Z^3*y-x^3*y+3*x^2*z)^2)

So are there 3 real-valued solutions on a common domain? Probably not.
Notice that the 3 solutions are not real or imaginary at the same xyz values.
Example:
 

sols:=allvalues(sol);
evalf(eval([sols],{x=1,y=2,z=3}));
simplify(fnormal(%));
evalf(eval([sols],{x=1,y=-2,z=3}));

limit(subs(theta2=theta1,soln[3]),theta1=0); #Returns unevaluated
limit(subs(theta2=theta1,soln[3]),theta1=0) assuming positive; #Undefined
limit(subs(theta2=theta1,soln[3]),theta1=0,right) assuming positive; #infinity
limit(subs(theta2=theta1,soln[3]),theta1=0,left) assuming positive; # -infinity

I think the responses are quite appropriate.

As Markiyan has already pointed out the two solutions are both correct. Yours is simpler though.
 

restart;
deq := diff(theta(x), x) = -sin(theta(x));                         
sol := dsolve({deq, theta(0) = -(1/2)*Pi});
sol:=simplify(expand(sol)) ;
solp:=simplify(sol) assuming x>0;
soln:=simplify(sol) assuming x<0;
sol2:= theta(x)=-2*arctan(exp(-x));
simplify(rhs(solp)-rhs(sol2)) assuming x>0;
simplify(rhs(soln)-rhs(sol2)) assuming x<0;

My favorite method for extracting solutions from solve (or dsolve) is as follows:

P := solve({3*y+x, 2*x+3+y});
X,Y:=op(subs(P, [x,y] )); 

Then you can use X and Y as you please, as in X+Y.
I don't ever assign to unknowns: notice the use of upper case X and Y.
You can if you want to though, assign to x and y. This can be done by just doing
 

assign(P);

Then try
x;   y;