1. You have 30 data values for each of S, P1, and P2_P2e. What are the corresponding times?
2. You say that "the initial point" is not known. By that I guess you mean the option initialpoint to LSSolve. By initialpoint in this sense is only meant a start guess just like in Newton's method you need a start guess for the root finding to get started. If not provided to LSSolve it will try to find one itself.
But you write that k1=0.000438, k4=0.0385, keq=2.7385 and T1_p1=36.8 are what you hope to get. Thus you have a fine initialpoint right there. 
3. Are the parameters beta and Tr in K:=t->cos((1/180)*beta*Pi)^(t/Tr) known? If so what are they?
4. It is somewhat strange to me that K as defined above should be just multiplied onto the ode solutions S,P1,P2,P2_e. So K doesn't enter into the ode system but is just applied afterwards? Any explanation for that?
 

@AlexSss If you are used to programming, then you might like to change the default input input from 2D Math Input to Maple Notation (aka 1D Math Input). At least you should try it.
If you do then space will not be interpreted as multiplication and you can write
with  (LinearAlgebra):
with as many spaces as desired.

To change go to  Tools/Options/Display/Input Display/  and change to Maple Notation.
Keep the Output display as it is.
While still in Tools/Options/ go to Interface/Default format for new worksheets. Change from Document to Worksheet.
Now click the button Apply Globally.

Changes take affect when you open a new worksheet.
These changes can be undone by the analogous procedure at any time.
I never use anything but Maple Notation and Worksheet mode.

@Harry Garst I should have mentioned that since the eigenvalue 0 has geometric multiplicity 4, the corresponding 4 eigenvectors seen in the last four columns of V might as well have been any other 4 linearly independent vectors spanning the same 4-dimensional space. They will in any case be orthogonal to the eigenvector representing the nonzero eigenvalue, but they need not be (and are in fact as we see in the symbolic case) not mutually orthogonal (and so not orthonormal). They could be made so by the Gram-Schmidt orthogonalization procedure.
## But using Gram-Schmidt in the symbolic case without first evaluating the symbolic results at the value of your parameters will create huge symbolic results even if you just try the last 4 columns
GramSchmidt([seq(V[..,i],i=2..5)],normalized);
But if you evaluate V at parameters then GramSchmidt is no problem:

 

params:={lambda[1,1] = 0.9,lambda[2,1] =0.8,lambda[3,1] =0.7,lambda[4,1] = 0.85,lambda[5, 1] = .7, theta[1,1] = 0.25,theta[2,2] = 0.21,theta[3,3]= 0.20,theta[4,4] = 0.15,theta[5,5] = 0.35};
V1:=eval(V,params);
GS:=GramSchmidt([seq(V1[..,i],i=1..5)],normalized);
V1n:=Matrix(GS);
fnormal(V1n^%T.V1n); #The identity matrix
OM1:=map(eval,Omega,params);
simplify(fnormal(V1n^%T.OM1.V1n));

@student_md You can save yourselt trouble of this sort by using a slightly different notation for (in this case) X[1] and w[1] in the actual computations. But that may be awkward. Delaying evaluation by using '  ' only does just that. Even in this case because the axis labels are kept in the plot structure and if that is evaluated again, then the assignment to X[1] (or w[1]) shows up instead.
Below I show a device that goes the other way: Changing the name of the label in an invisible way.
Replace w[1] by  ` w`[1] . That is backquote, space, w, backquote.

restart;
w[1]:=67: #Anything will do
plot1:=plot(x^(1/2),x=0..1,labels=[c,typeset('w[1]')]);
plot1; #67 becomes label now
plot1:=plot(x^2,x=0..1,labels=[c,typeset(` w`[1])]); #The little trick
plot2:=plot(x^(1/2),x=0..1,color=blue); #  If the labels are the same no need for them here.
plots:-display(plot1,plot2);

Finally I should mention that by using backquotes you can turn "anything" into a name, as in this silly example, where the name is `56`.
`56`:=123;
 

@sideshow To see what is going in that expression, try the same thing with an unassigned function F.
 

T:=add((D@@i)(F)(0)*m^i/i!,i=0..2);

D is the operator of differentiation of functions (procedures). So D(sin) returns cos.
D(sin)(Pi/2) therefore returns 0 since cos(Pi/2)=0.
(D@@2)(F) is the second derivative of F, so (D@@2)(sin) = -sin, and so (D@@2)(sin)(Pi/2) = -1.
For convenience (D@@0)(F) is simply F (or if you will, the 0th derivative of F).
That is used above in the add command where I have i running from 0 to 2.
## An ugly looking alternative (in my view) is this:
 

Td:=add(eval(diff(F(m),[m$i]),m=0)*m^i/i!,i=0..2);

For concrete functions like the function f defined in my answer above, T and Td return the same:
 

f:=m->tanh(sqrt(q)*z/y+x*m/y);
eval(T,F=f);
eval(Td,F=f);

You may have a look at the SignalProcessing package:
?SignalProcessing
I don't know anything about that, but maybe you will find what you need there.

Try removing expand (only those letters).  Assign the resulting expression to u, say.
Then do:
 

eval(u,{a=1,h=2,o=3,v=4});
eval(expand(u),{a=1,h=2,o=3,v=4});
w:=-2*a^23*h^40*o*v^2+2*a^23*h^40*v^3;
eval(w,{a=1,h=2,o=3,v=4});

You should find that the first 2 results are the same, but the result you report as the Maple 17 GUI and Mathematica result (w) is WRONG.
Also notice that I cannot confirm the difference you report between Maple 17 GUI and command line.
I use Maple 17.02.

@gkokovidis I tried the GUI and the command line version in Maple 17.02.
I got the same result as you report in both cases.
The same is true of Maple 2016.2.

@rath1 Yes, I expected that answer.
I my opinion it was a very bad mistake, when 2D-math input was introduced, to make a space mean multiplication. 
In particular in Maple, where parentheses () are used for two purposes:
(1) for function application as in f(x)
(2) for grouping of terms as in a*(b+c)
In Mathematica they can get away with it since function application uses square brackets:  f[x].

I'm not personally bothered by the problem since I use Maple notation (aka 1D-math) exclusively for input.
Let me add that because I use a Danish keyboard, I'm pleased with parentheses in function application, since the square brackets need Alt Gr + [ .
But then again square brackets are used extensively in Maple for other purposes.

@esmerindobernardes I'm just saying that your commands as I gave them return what seems to me to be a perfectly good answer. If I had thought I was really helping, I would have called this an answer!

For your information here are the details of my installation:

kernelopts(version);
   Maple 18.02, X86 64 WINDOWS, Oct 20 2014, Build ID 991181
interface(version);
 Standard Worksheet Interface, Maple 18.02, Windows 8, October 20 2014 Build ID 991181

My executed worksheet is here: MaplePrimes2017-02-24_Maple18.mw

I suppose that you won't be satisfied with using Maple Notation as output?
You can change to that in Tools/Options/Display/Output Display.
 

I tried your code in my Maple 18.02.
 

restart;
Hso := Matrix(8, {(1, 4) = -x, (1, 6) = I*x, (2, 3) = x, (2, 5) = I*x, (3, 2) = x, (3, 5) = -I*z, (3, 8) = y, (4, 1) = -x, (4, 6) = I*z, (4, 7) = -y, (5, 2) = -I*x, (5, 3) = I*z, (5, 8) = -I*y, (6, 1) = -I*x, (6, 4) = -I*z, (6, 7) = -I*y, (7, 4) = -y, (7, 6) = I*y, (8, 3) = y, (8, 5) = I*y});
LinearAlgebra[Eigenvectors](Hso);

No problem. Large output.

@shesia A good book is the Programming Guide, which is available as a series of worksheets and as a pdf file.
In Maple itself you can access these worksheets by doing:
? Programming Guide

You can also download a pdf version from
http://www.maplesoft.com/documentation_center/

@shesia OK, I tried with the much larger set of data.
NOTE: Are these data really measured data? With that much accuracy! Or are they constructed by a teacher?
I revised Q slightly to make it:
(1) catch possible errors and in that case return the large list of residuals [10^6$(3*nops(T1))].
(2) print the parameter points for which Q is computed to show that something is going on.
 

Q:=proc(T1_p1, k1, keq, k4) option remember; local P1v,P2v,P2e_v,Sv,resid;
   res(parameters=[T1_p1, k1, keq, k4]);
   printf("Trying %a\n",[_passed]);
   try
     P1v:=P1fu~(T1);
     P2v:=P2fu~(T1);
     P2e_v:=P2e_fu~(T1);
     Sv:=Sfu~(T1);
     resid:=[P1v-P1d1,P2v+P2e_v-P2_P2ed1,Sv-Sd1];
     return [seq(seq(resid[i][j],i=1..3),j=1..nops(T1))]
   catch:
     return [10^6$(3*nops(T1))]
   end try
end proc: 

Then with q defined as

q:=[seq(subs(_nn=n,(proc(T1_p1, k1, keq, k4) Q(_passed)[_nn] end proc)),n=1..3*nops(T1))]: 

and the command

solLS:=Optimization:-LSSolve(q,initialpoint=[15,0.02,4,4]);

I got (after a while)

The plots don't look as good as before: That is not surprising if the data are real data (i.e. not constructed by a teacher as I remarked about above).
After all, isn't the differential equations model you have only a crude model?
A thing to keep in mind is that LSSolve doesn't necessarily find the global minimum!
For the same plot as I reported for the small data set I got:

@shesia The very simple reason is that you have 51 t-values, but only 24 values in the other lists:
nops~([T1,Sd1,P1d1,P2_P2ed1]);
                                                              [51, 24, 24, 24]