Copying your ode I find that you have a sequence of length 2 on the left (and a zero on the right).
Your ode is:

ode:=((0.9181604810e-6*cos(phi(t))^3+0.4213239281e-4*cos(phi(t)))*(diff(diff(phi(t), t), t))-0.9181604810e-6*(diff(phi(t), t))^2*cos(phi(t))^2*sin(phi(t))+0.2166831514e-7*(diff(phi(t), t))^2*sin(phi(t))^3*cos(phi(t))^2/(-0.2359970352e-1*cos(phi(t))^2+1)-0.5976753198e-5*(-.1536219500*(diff(diff(phi(t), t), t))*sin(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)-.1536219500*(diff(phi(t), t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.3625432474e-2*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(3/2))*sin(phi(t))*(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)*cos(phi(t))+0.205427606e-5*(-.1536219500*(diff(diff(phi(t), t), t))*sin(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)-.1536219500*(diff(phi(t), t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.3625432474e-2*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(3/2))*cos(phi(t))+0.5376768250e-3*(diff(phi(t), t))*cos(phi(t)), 

0.1836916543e-3*cos(phi(t))^2*(diff(diff(phi(t), t), t))-0.1836916543e-3*sin(phi(t))*cos(phi(t))*(diff(phi(t), t))^2-0.2821907012e-4*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.7783642454e-2*(-.1536219500*(diff(diff(phi(t), t), t))*sin(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)-.1536219500*(diff(phi(t), t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.3625432474e-2*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(3/2))*(-0.2359970352e-1*cos(phi(t))^2+1)+0.1210411733e-2*(diff(diff(phi(t), t), t))*sin(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)+0.1210411733e-2*(diff(phi(t), t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(1/2)-0.2856535802e-4*(diff(phi(t), t))^2*sin(phi(t))^2*cos(phi(t))/(-0.2359970352e-1*cos(phi(t))^2+1)^(3/2)+1.553675255*(50*exp(-2.601456816*(phi(t)-4*Pi*trunc((1/4)*(phi(t)+Pi)/Pi)-.36)^2)+2)*Pi*cos(phi(t))) = 0;

Surely that comma ougth to be something else or you only want one of  the two?
The error comes in at the definition of Bewegungsgl:
nops([lhs(Bewegungsgl)]); # answer 2
In fact already here:
nops(GleichungKräfte[2]);

You can do as follows:

restart;
eq:=diff(x(t), t$2) = .8/(x(t)^3*exp(0.02*int(1/x(s), s = 0 .. t-tau)))-1/x(t)^2;
ics := x(0) = 1, D(x)(0) = 0;
eq_y:=y(t)=int(1/x(s), s = 0 .. t-tau);
ic_y:=eval(eq_y,{t=0,x=1});
sys:={diff(eq_y,t),subs(rhs(eq_y)=lhs(eq_y),eq)};
res:=dsolve(eval(sys union {ics,ic_y},tau=0.1),numeric);
plots:-odeplot(res,[t,x(t)],0..30);

After giving your system the name sys I did:
 

g:=indets(sys)=~1;
fsolve(sys,g);

and quickly got
{c[0][0] = 12.89605073, c[0][1] = -6.470667865, c[0][2] = 1.752111454, c[0][3] = -.1774943190, c[1][0] = 9.737598780, c[1][1] = -2.590154874, c[1][2] = .3402465368, c[1][3] = -0.1700267263e-1, c[2][0] = 7.032179265, c[2][1] = -1.179016309, c[2][2] = 0.9696822774e-1, c[2][3] = -0.3121342845e-2}
which agrees well with your second real solution.
The title of your question says that fsolve gives an error with a range. It didn't in your worksheet; it just returned unevaluated, i.e. it didn't succeed in finding a solution.
Often I find that using a guess instead of a range (or ranges) is surer of success.

Even the version in the help page for semitorus in Maple 2016.2 looks weird.
Using capped=false gives you a nice plot, but with no cap.
Here is what I got in Maple 2015.2 using your command as given, i.e.

plots:-display(plottools:-semitorus([0, 0, 0], 0 .. Pi, 1, 2), lightmodel = light4, orientation = [-140, 60], scaling = constrained, style = patchnogrid);

If you turn the plot (in Maple, not here) you will see that it looks like a donut cut in half horizontally and with a piece of paper at the bottom. The "paper" is due to the default capped=true.

I get the same thing in Maple 18 and 15, while in Maple 12 I get something quite different:

Obviously here the angle means something different than in the other versions.

To get something like this Maple 12 version in Maple 2016.2 you could try tubeplot (in the plots package. Here I take 3/4 of the donut cut vertically, but not capped:

tubeplot([cos(t),sin(t),0],t=0..3*Pi/2,radius=1/2,numpoints=60,tubepoints=20,scaling=constrained);

I just checked in Maple 8. It works as Maple 12.
I wonder if there was a documented change in Maple 13 or 14 about the meaning of the angular range (the second argument) in semitorus.
My suspicion is that this is a bug introduced in Maple 13, 14, or 15.
I notice that the help pages are the same.
In Maple 2016 you find
"The semitorus command creates a three-dimensional plot data object, which when displayed is a section of the torus centered at [x,y,z], whose meridian has radius r and distance from the center of a meridian to the center of the semitorus is R.  The semitorus starts at radian angle a and ends at radian angle b."
In Maple 12 you find:
"The semitorus command creates a three-dimensional plot data object, which when displayed is a section of the torus centered at [x,y,z], whose meridian has radius r and distance from the center of a meridian to the center of the semitorus is R. The semitorus starts at radian angle a and ends at radian angle b."

That looks the same to me.
I will submit an SCR.

In a recent version of Maple (I tried Maple 12, 18, 2016 as a sample) try:
'rand()/rand()'; # You get one
rand()/rand(); #Highly unlikely that you get one
If you try the same two commands in Maple 8 you get 1 in both cases.
## You could try the difference also:
##In Maple 8 both of these give 0:
'rand()-rand()';
rand()-rand();
## whereas in Maple 12, 18, and 2016 you only get 0 for the first.
## If you look at the code for rand in recent versions (including 12) you will see that
RandomTools:-MersenneTwister:-GenerateInteger is used.
If you look at the code for rand in Maple 8, you will find it totally different. (use showstat(rand)).
###
While the rand procedures are surely different the reason for the different behavior in Maple 8 and Maple 2016 (say) is not that the procedures are written differently.
I tried copying the procedure from Maple 8 into Maple 2016 and gave it the name rand8.
Result:
rand8()/rand8(); #doesn't return 1 in Maple 2016.
###############
Now try this:

_n:=1;
p:=proc() global _n; _n:=_n+1; _n end proc;
p()/p(); #returns 2/3 in Maple 2016, but 1 in Maple 8
'p()/p()'; #returns 1 in both versions

Thus the change from Maple 8 to more recent versions has to do with the handling of automatic simplification.
## It should be added that p()+p() returns the same as 2*p() would, i.e. a sum like a+a is automatically simplified to 2*a.

Numerical solution has been available for quite a while.
Example with L=2:
 

restart;
ode:=diff(y(x),x$2)+lambda*y(x)=0;
bc:=y(0)=0,y(L)=0;
# Taking as an example L=2:
res:=dsolve(eval({ode,bc,D(y)(0)=1},L=2),numeric);
res(0); # You see lambda
plots:-odeplot(res,[x,y(x)],0..2);

The result of the command Sol2 := pdsolve({IBC2s, PDE2});  is NULL, i.e. Maple can't do it.
I'm using Maple 2016.2.
Try
evalb(Sol2=NULL); # true
So you are feeding build nothing, as in
build();
###
Try without IBC2s:

Sol3 := pdsolve(PDE2);
build(Sol3);

I tried your code in Maple 2016.2, Maple 2015.2, Maple 18, and Maple15:
 

restart;
with(Statistics):
X := RandomVariable(StudentT(nu));
PDF(X,0.5);

and got:

.5641895835*GAMMA(.5000000000*nu+.5000000000)/(sqrt(nu)*GAMMA(.5000000000*nu)*(1.+.25/nu)^(.5000000000*nu+.5000000000))

So I have no idea what is going on.

You will surely get some better answer, but here is a very crude way:
Insert a computation that takes a very long time at the place desired.
Interrupt the computation by using the button provided for that.
Example:
 

restart;
a:=8;
b:=89;
##
for i to 10^16 do evalf(sin(i)) end do:
##
c:=56;
d:=324;

Well, I got an answer, not great, but the supposedly good values given in your question don't seem to better at all.
 

restart;
ode_sub := diff(S(t), t) = -k1*S(t)-S(t)/T1_s;
ode_P1 := diff(P1(t), t) = 2*k1*S(t)-k2*(P1(t)-P2(t)/keq)-P1(t)/T1_p1;
#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;
ode_system := ode_sub, ode_P1, ode_P2, ode_P2e;
s0 := 96304.74567; k2 := 10^5; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1;
init := S(0) = s0, P1(0) = 0, P2(0) = 0, P2_e(0) = 0;
ode_system;
T := [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];
####
data_s := [96304.74567, 77385.03700, 62621.83067, 51239.94333, 42663.82367, 35084.74100, 28480.28367, 23066.01467, 18774.73700, 15179.13700, 12278.50767, 9937.652000, 8046.848333, 6521.242000, 5287.811667, 4277.779000, 3466.518333, 2835.467000, 2297.796333, 1861.249667, 1529.654000, 1235.353000, 999.6626667, 826.2343333, 667.9480000, 559.9230000, 449.2790000, 376.4860000, 289.1203333, 236.1483333];
####
data_p1 := [0.86e-1, 3.904, 26.975, 31.719, 41.067, 46.779, 52.115, 43.101, 44.344, 41.094, 36.523, 27.742, 26.543, 28.062, 22.178, 21.303, 14.951, 17.871, 11.422, 12.051, 9.232, 6.817, 6.1, .717, 1.215, 6.146, .772, .375, 2.595, .518];
####
data_p2_p2e := [-3.024, 22.238, 61.731, 103.816, 132.695, 159.069, 167.302, 160.188, 158.398, 152.943, 146.745, 135.22, 132.145, 120.413, 107.864, 95.339, 90.775, 81.828, 71.065, 70.475, 62.872, 49.955, 40.858, 42.938, 41.311, 35.583, 31.573, 29.841, 29.558, 21.762];
#### The data are vastly different in size, so it might be good to put weights on the residuals.
#### I have done that in the uncommented version of resid in QM
max(data_s);
max(data_p1);
max(data_p2_p2e);
####
beta:=10: Tr:=2:
K:=unapply(evalf(cos((1/180)*beta*Pi))^(t/Tr),t);
KT:=<K~(T)>;
####
resM:=(T1_p1xx, k1xx, keqxx, k4xx)->
  dsolve(eval({ode_system,init},{T1_p1=T1_p1xx, k1=k1xx, keq=keqxx, k4=k4xx}),numeric,_rest);
#### Dividing the data elementwise by the elements in KT:
P1d:=<data_p1>/~KT;
P2_P2ed:=<data_p2_p2e>/~KT;
Sd:=<data_s>/~KT;
###############################################################################
QM:=proc(T1_p1, k1, keq, k4) option remember; local res,A,P1v,P2v,P2e_v,Sv,resid;
    userinfo(1,procname,printf("Trying %a\n",[_passed]));
    try
     res:=resM(T1_p1, k1, keq, k4,maxfun=10^8,output=Array(T),stiff=true); #stiff used
     A:=res[2,1];
     ## 
     if A[-1,1]=`?` then error "cannot evaluate the solution" end if;
     catch "cannot evaluate the solution":
       return [10^6$(3*nops(T))];
     end try;
     P1v:=A[..,2];
     P2v:=A[..,3];
     P2e_v:=A[..,4];
     Sv:=A[..,5];
     #resid:=[P1v-P1d,P2v+P2e_v-P2_P2ed,Sv-Sd];# No weights
     resid:=[(P1v-P1d)/max(P1d),(P2v+P2e_v-P2_P2ed)/max(P2_P2ed),(Sv-Sd)/max(Sd)];#Weight put on all
     [seq(seq(resid[i][j],i=1..3),j=1..nops(T))]
     end proc:
q:=[seq(subs(_nn=n,(proc(T1_p1, k1, keq, k4) QM(_passed)[_nn] end proc)),n=1..3*nops(T))]: 
#################################################################################
## The wrapping of Codetools:-Usage around the LSSolve command is done for measuring time etc.
##
forget(QM):
sol:=CodeTools:-Usage(Optimization:-LSSolve(q,initialpoint=[0.03, 0.02, .5, .5],output=solutionmodule));
sol:-Results();
sol:-Settings();
Number of calls to QM:
nops([indices(op(4,eval(QM)))]);
params:=convert(sol:-Results(solutionpoint),list);
RES:=resM(op(params),stiff=true);
pS:=plot(T,data_s,style=point,symbolsize=20,symbol=solidcircle):
pP1:=plot(T,data_p1,style=point,symbolsize=20,symbol=solidcircle):
pP22e:=plot(T,data_p2_p2e,style=point,symbolsize=20,symbol=solidcircle):
plots:-display(pS,plots:-odeplot(RES,[t,S(t)],0..T[-1]),labels=[t,S]);
plots:-display(pP1,plots:-odeplot(RES,[t,P1(t)],0..T[-1]),labels=[t,P1],thickness=3); 
plots:-display(pP22e,plots:-odeplot(RES,[t,P2(t)+P2_e(t)],0..T[-1]),labels=[t,P2+P2_e]);
The values given as "the solution" don't fit very well:
RESsol:=resM(36.8,0.000438,2.7385,0.0385,stiff=true);
plots:-display(pS,plots:-odeplot(RESsol,[t,S(t)],0..T[-1]),labels=[t,S]);
plots:-display(pP1,plots:-odeplot(RESsol,[t,P1(t)],0..T[-1]),labels=[t,P1],thickness=3); 
plots:-display(pP22e,plots:-odeplot(RESsol,[t,P2(t)+P2_e(t)],0..T[-1]),labels=[t,P2+P2_e])

As the setting is above the result is
params := [0.0501239039106269, 0.0115045966178040, 0.663492002313513, 0.434185918627229]
Remember that the order is [T1_p1, k1, keq, k4].
The graphs are not great. My only comment to that is one I have given before:
If the data are real data, then it is unlikely that your model is all that good. After all it is linear, so without knowing anything about the origin of the model, my guess is that it is used because it is the simplest one could come up with. Never a bad start though.

If weights are not put on the residuals the graph for S looks almost perfect, but that is no miracle since data are so big.

First of all you can only work on a finite interval. Thus you must replace infinity by some (in the context) large number.
I have used subs(infinity=4,{bc}) below, but you can just remove infinity and write 4 instead.
Secondly, the second derivative in bc should be written (D@D)(f)(0) or alternatively as (D@@2)(f)(0). The latter version is easier to use for higher derivatives.
Thirdly, cosine is spelled cos, not Cos.
Fourthly, use a midpoint method.
After these corrections the following works. Notice that I use infinity=4.

res:=dsolve({eqn1,eqn2,eqn3} union subs(infinity=4,{bc}),numeric,method=bvp[midrich]);
plots:-odeplot(res,[eta,f(eta)]);

Try this:
 

f:=x->piecewise(0<x and x<0.5,x,0.5<x and x<1,1-x, 0);
g:=unapply(convert(f(x),Heaviside),x);
# Then compare:
diff(f(x),x);
diff(g(x),x);

Just use evalc as in
v := evalc(Vector(4, [Re(V[1, 1]), Im(V[1, 2]), Re(V[2, 1]), Im(V[1, 1])]));

evalc assumes itself that variables are real.
No need for making assumptions yourself.
You had a syntax error: When I copy the assume line I get
`assume `*{psi(t), 'real'}
which is utter nonsense:
(1) Curly braces are used for defining sets only.
(2) The appearance of the multiplication sign is due to the silly convention in 2D math input (the unfortunate default) that a space means multiplication.
But the good news: You don't need assume here.

You could normalize the rows like this:

M := Matrix([[1,2],[3,4]]);
<seq(M[i]/LinearAlgebra:-Norm(M[i],2),i=1..2)>;

Ij you have an initial value problem you can set abserr and relerr.
?dsolve,numeric,ivp  and the reference in that page to Error Control for Numerical Solution of IVPs in Maple
If you have a boundary value problem you can set abserr.
?dsolve,numeric,bvp