About computing integral integrate(exp(-x^k),x,0,1)

About computing integral   integrate(exp(-x^k),x,0,1)

(%i1) S(k):='integrate(exp(-x^k),x,0,1)$

 Note: output from workshet is removed. Plese copy to wxMaxima and execute.

(%i2) assume(x>0);
 Correct results can be obtained in the following ways:

 a) using formula
integrate(f,x,0,1)=integrate(f,x,0,inf)-integrate(f,x,1,inf):
(%i3) for k thru 10 do
(integrate(exp(-x^k),x,0,inf)-integrate(exp(-x^k),x,1,inf),
print(S(k)=%%,"~~",float(%%)));

 b) using Newton-Leibniz formula  integrate(f,x,a,b) = F( b)-F(a) :
(%i4) for k thru 10 do
(integrate(exp(-x^k),x),
subst(x=1,%%)-subst(x=0,%%),
float(%%),
print(%%));

 c) using change of variable y=-x :
(%i5) for k thru 10 do
(changevar(S(k),y=-x,y,x),ev(%%,nouns),
print(S(k)=%%,"~~",float(%%)));

 d) using change of variable y=1/x :
(%i6) for k thru 10 do
(changevar(S(k),y=1/x,y,x),ev(%%,nouns),
print(S(k)=expand(%%),"~~",float(%%)));

 e)  using change of variable x=y^(2/k) :
(%i7) for k thru 10 do
(changevar(S(k),x=y^(2/k),y,x),ev(%%,nouns),
print(S(k)=expand(%%)));

 f)  using change of variable x=sin(y):
(%i8) for k thru 10 do
(changevar(S(k),x=sin(y),y,x),ev(%%,nouns),
print(S(k)=expand(%%)));

 g)  using change of variable x=cos(y):
(%i9) for k thru 10 do
(changevar(S(k),x=cos(y),y,x),ev(%%,nouns),
print(S(k)=expand(%%)));

  But, Maxima compute wrong:
(%i10) for k thru 10 do
(integrate(exp(-x^k),x,0,1),
print(S(k)=%%));

 I think , that Maxima  wrong compute like this:
(%i11) for k thru 10 do
(a:changevar(S(k),exp(-x^k)=y,y,x),
changevar(%%,-log(y)=z,z,y),
b:factor(%%),
ev(%%,nouns),
print(((S(k)=a)=b)=%%));

 Aleksas D