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;

Starting from the odes and leaving out the data corresponding to initial conditions (they won't help)  we can do:
 

restart;
ode1 := diff(SS(t), t) = -k*SS(t)-SS(t)/B;
ode2 := diff(PP(t), t) = k*SS(t)-PP(t)/T1p;

init:=PP(0)=0,SS(0)=100000;
dsolve({ode1,ode2,init});
sol:=combine(expand(%));
PS:=subs(sol,[PP(t),SS(t)]);
T:=<2,4,6,8>:
S:=<86089.76983, 74114.4849, 63804.98946, 54929.56828>:
P:=<8932.499583, 15965.05975, 21410.35044, 25533.98716>:
###
Pfu,Sfu:=op(map(unapply,PS,t));
RP:=convert(P-Pfu~(T),list);
RS:=convert(S-Sfu~(T),list);
RPS:=[op(RP),op(RS)];
###
res:=Optimization:-LSSolve(RPS,k=0..1,B=0..1000,T1p=0..10);
SOL:=eval(sol,res[2]);
Psol,Ssol:=op(subs(SOL,[PP(t),SS(t)]));
plot(T,P,style=point,symbolsize=20); p1:=%:
plot(T,S,style=point,symbolsize=20); p2:=%:
plot(Psol,t=0..8,color=blue); p1s:=%:
plot(Ssol,t=0..8,color=blue); p2s:=%:
plots:-display(p1,p1s);
plots:-display(p2,p2s);

You have some syntax problems. I have corrected them below. Besides, X has 6 elements while Y has 5. I removed 0 from X.
By the assignment model:= 1.*10^5*exp(-t*k-*t/B)    I have assumed that you meant  model:= 1.*10^5*exp(-t*k-t/B).

Then I observe that exp(-t*k-t/B) has really only one parameter, namely  a=k+1/B.
Thus you can only expect unique results for the quantity k+1/B not for k and B individually.
## Maybe you wanted

model:= 1.*10^5*exp(-t*k-B/t);

which actually gives a good fit. ?

restart;
model:= 1.*10^5*exp(-t*k-t/B); #corrected
X:=<2,4,6,8,10>; #corrected
Y := <100000, 86089.76983, 74114.4849, 63804.98946, 54929.56828>;
res:=Statistics:-Fit(model,X,Y,t);
Statistics:-Fit(model,X,Y,t,output=[parametervalues]);
eval(model,op(%));

There is agreement between the last output and the first output from Fit.
### Is your mdel actually what I corrected it to be or what is it?
Plotting:
 

plot(X,Y,style=point, symbolsize=20); p1:=%:
plot(res,t=2..10,color=blue); p2:=%:
plots:-display(p1,p2);

I made a procedure way back when:
 

`convert/sn`:=proc(f) local SN;
SN:=t->if ilog10(t)=0 then t else 10^(-ilog10(t))*t*10^``(ilog10(t)) end if;
subsindets(f,float,SN)
end proc:
convert(evalf(1000*Pi),sn);
expand(%);

Express the right hand side of the last line with the help of piecewise:
p:=piecewise(m>=0, cos(m*theta), -sin(m*theta));