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]

I noticed that when I copied your code s0 was set at 100000 which doesn't agree very well with the data.
I notice now that s0 is set at 10000, which is fine. I suppose you edited that in the meantime?
I shall edit my answer below to that effect.

In Maple 2016 the error message for your example (with D[1] as you have it) is more informative:

Error, (in PDEtools:-Solve) unable to handle functions of not names as in the case of the system passed, containing [_F2(-x)]

To see why there is talk about _F2(-x) try

pdsolve(PDE);
eval(%,t=0);

In Maple 12 (I don't have 13) the error message is

Error, (in pdsolve/sys) unable to handle PDE problem subject to boundary conditions {u(x, 0) = cos((1/2)*x), (D[1](u))(x, 0) = 0}

You don't need the package PDEtools, but if you do, the syntax is
with( PDEtools );

@Thomas Dean Yes, if D[2] is meant (i.e. the first derivative with respect to the second variable), then there is no problem:

pdsolve({PDE, u(x, 0) = cos((1/2)*x), D[2](u)(x, 0) = 0});

The warning just means that solve cannot find any (or all) solutions, but on the other hand cannot claim that there aren't any either.

If you want somebody to try your code, you need to upload a worksheet.

Your "code" doesn't tell us what npar, de, nde, Ctrl (and more?) are, so please upload a worksheet or, alternatively, give us the full code as text.

@Axel Vogt As for your new polynomial isn't Digits=15 way too low?
This works perfectly:
 

restart;
N:=40:
Digits:=100:
F:=expand(mul(I*x-k,k=1..N)):
f:=evalf(F):
S:=[fsolve(f,complex)];