The presence of the imaginary part after the use of evalf is due to rounding. Try evalf[20](...) .

restart;
s1 := (1/(2*Pi)+Q*cos(k*phi))^2/(1+(1/(2*Pi)+Q*cos(k*phi))^2);

Here is my version built upon yours.
Q:=15;
#The integral from 0 to 2*Pi is what you want and that is independent of the integer k, so:

k:=1;
#The integral over 0..2*Pi is twice the integral over 0..Pi
A:=Int(s1,phi=0..Pi);
evalf(A);
#Changing variables.
IntegrationTools:-Change(A,phi=2*arctan(u),[u]);
evalindets(%,specfunc(anything,cos),simplify@expand);
normal(%);
s2:=IntegrationTools:-GetIntegrand(%);
t:=singular(s2,u);
t1a:=map(rhs@op,[t]);
t1:=select(x->evalb(evalf(evalc(Im(x))>0)),t1a);
Res:=add(residue(s2,u=p),p=t1):
simplify(2*Pi*I*Res);
evalf[20](%);

If on the surface of it you see common factors you could force factorization by a procedure like the following.

sfactor:=proc(u) local i,F,S,p;
if not hastype(u,`+`) then return u end if;
if not type(u,`+`) then return evalindets(u,`+`,procname) end if;
for i to nops(u) do
  if not type(op(i,u),`*`) then
    F[i]:={op(i,u)}
  else
    F[i]:={op(op(i,u))}
  end if
end do;
S:=`intersect`(seq(F[i],i=1..nops(u)));
if S={} then
   u
else
   p:=`*`(op(S));
   p*map(`/`,u,p)
end if
end proc;

#Test on your example.
u:=(Pi^2*(c+(1+Pi)/Pi)^2+Pi^2)/((c-(-1+Pi)/Pi)^2*Pi^2+Pi^2);
T:=op(1,u);
sfactor(T);
sfactor(u);

What takes time is evaluating the sum. When k = 0 that time is wasted, since the result in your particular example obviously is zero.

In your example you could have proceeded in one of the following two ways.

#Method 1: unevaluation quotes around s.
restart;
f:=(k,y)->k*'s'(y);
s:=y->add(sin(1/y^i),i=1..1000000);
f(0,y);
f(1,y);
%;
#Method 2: using the inert form for add, %add.
restart;
f:=(k,y)->k*s(y);
s:=y->%add(sin(1/y^i),i=1..1000000);
f(0,y);
f(1,y);
value(%);

With no details about the more general cases your mention it is impossible to say if these ideas are useful.

In the help page for SingularValues you find:

"The singular values S together with the left singular Vectors as columns of Matrix U and the right singular Vectors as columns of Matrix Vt satisfy
                       (U . SS) . Vt = A
. Matrix SS has the entries of S along its main diagonal."

Example:
restart;
with(LinearAlgebra):
A := RandomMatrix(5,3,datatype=float);
S,U,Vt:=SingularValues(A,output=['S','U','Vt']);
S1:=DiagonalMatrix(S[1..3],5,3);
U.S1.Vt;
Equal(evalf(%),A);

If you are not already aware of dualaxisplot from the plots package, you may want to explore the possibilities it offers together with options for label directions and gridlines.

You can try removing the initial condition:

remove(has,edeq,_C[0]);
dsolve(%);

You will see (this is Maple 12 at the moment) that Maple cannot find an explicit expression for the general solution, so it is no wonder that solve must return NULL when attempting to solve the equation involved in satisfying the initial condition.

I assume that you meant

A:= subs(Fsubs,parvals,Theta0/C=(exp(Pec_i*F*(1-S0))-1)*(H+1/(1-exp(-C/k*Pec_i*F*S0))));

B is not an equation and seems to be missing a parenthesis. If Maple expects an equation and don't get one, often it assumes you meant equal to zero. Probably you didn't mean that.

With the modification of A above you have a connection between R3  and S0. You can explore that by doing

plots:-implicitplot(A,R3=-3..3,S0=-2..2,gridrefine=2);

Suppose you are interested in the case S0 >0, then you could define numerically S0 as a function of R3 like this:

f:=r->fsolve(eval(A,R3=r),S0=1);

#Test and plot
f(0.12345);
plot(f,0..3);

As long as you have an initial value problem you shouldn't have much of a problem. With the conditions at tf you have an initial value problem.

I have access to Maple 12 right now so will make use of map instead of the elementwise operation ~ which is available in Maple 13 and later.

Using MatrixExponential from the LinearAlgebra package:

with(LinearAlgebra);
A:=Matrix([[-1/R/C,1/C,0,0],[-1/L,0,0,0],[-2,0,1/R/C,1/L],[0,0,-1/C,0]]);
b:=<0,0,2*nf,0>;
Eigenvalues(A);
W:=unapply(MatrixExponential(A,t),t):
X:=W(t-tf).<x1f,x2f,p1f,p2f>+map(int,W(t-s).b,s=tf..t);

Using dsolve (and making use of A and b above.

xp:=<x1(t),x2(t),p1(t),p2(t)>;
S:=map(diff,xp,t)=A.xp+b;
sys:=convert(lhs(S)-rhs(S),set);
#dsolve without specifying the method gives a messy result:
res1:=dsolve(sys union {x1(tf)=x1f,x2(tf)=x2f,p1(tf)=p1f,p2(tf)=p2f}):
# method=laplace is quicker and neater.
res2:=dsolve(sys union {x1(tf)=x1f,x2(tf)=x2f,p1(tf)=p1f,p2(tf)=p2f},{x1(t),x2(t),p1(t),p2(t)},method=laplace);

In Maple 15 the following works.

restart;
sys:={diff(b1(t), t) = I*u1(t),
diff(b2(t), t) = I*u2(t)-3*b1(t)/(2*omega),
diff(b3(t), t) = I*u3(t),
diff(u1(t), t) = I*b1(t)+2*u2(t)/omega-gradQ1,
diff(u2(t),t)=I*b2(t)-1/(2*omega)*u1(t) - gradQ2,
diff(u3(t), t) = I*b3(t)+I*b3(t)-gradQ3};
gradQ2:= (3*K*u2(t)-2*u1(t)*omega)/((9/4)*K+omega^2+omega^2/kappa^2);
gradQ1 := 1.5 * K * gradQ2/omega;
gradQ3 := gradQ2/kappa;
K := (x + t);
inits:={b1(0) = exp(-(1/25)*ln(2)*x^2), b2(0) = (3/25*I/omega)*ln(2)*x*exp(-(1/25)*ln(2)*x^2), b3(0) = 0, u1(0) = 0, u2(0) = 0, u3(0) = 0};
res:=dsolve(sys union inits,numeric,parameters=[x,kappa,omega],range=0..16,output=listprocedure);
B1,B2,B3,U1,U2,U3:=op(subs(res,[b1(t),b2(t),b3(t),u1(t),u2(t),u3(t)]));
#Test
res(parameters=[1,.125,.125]);
B1(2);
#Now we make a procedure accepting values for the 3 parameters and for t.
#It has been made to return unevaluated if the parameters are not of type numeric.
#B1(t) itself will return unevaluated if t is not of type numeric.
Qb1:=proc(xx,k,w,t) if not type([xx,k,w],list(numeric)) then return 'procname(_passed)' end if;
   res(parameters=[xx,k,w]);
   B1(t)
end proc;
#Test
Qb1(x,.125,.125,t);
Qb1(1,.125,.125,t);
Qb1(1,.125,.125,2);
with(plots):
complexplot(Qb1(1,.125,.125,t),t=0..1);
animate(complexplot,[Qb1(x,.125,.125,t),t=0..0.5,thickness=3],x=-1..1,frames=10);

#When you assigned 'res' above you may have noticed that Maple didn't print anything to the screen.
This seems like a bug to me since it will normally do that. What exactly triggers this I don't know, but have noticed that the the following command, where I is  simply replaced by 1, prints the expected.

Res:=dsolve(eval(sys union inits,I=1),numeric,parameters=[x,kappa,omega],range=0..16,output=listprocedure);


#I also tried the whole code in Maple 12 and it didn't work, but that is not surprising since even the parametrized initial value problem (notice the presence of exp(a), not just a or e.g. a^2).

eq:=diff(y(t),t)=y(t):
res2:=dsolve({eq,y(0)=exp(a)},numeric,parameters=[a]);

#gives the error
Error, (in dsolve/numeric) invalid specification of initial conditions

Added. I have submitted a bug report (aka "Software change request"). The printing problem has nothing to do with using the parameters option, but seems to occur because of the presence of imaginary numbers in the ivp.

Here is a solution to the first system using the procedure option in dsolve/numeric
?dsolve,numeric,ivp

restart;
odes := diff(x(t), t, t)*sin(diff(y(t), t))+diff(y(t), t, t)^2 = -2*x(t)*y(t)*exp(-diff(x(t), t))+x(t)*diff(y(t), t)*diff(x(t), t, t), x(t)*diff(y(t), t, t)*diff(x(t), t, t)+cos(diff(y(t), t, t)) = 3*y(t)*diff(x(t), t)*exp(-x(t)):
ics := x(0) = 1, (D(x))(0) = 0, y(0) = 0, (D(y))(0) = 1:
#We rewrite the system as a first order system:
subs(diff(x(t),t,t)=yp2,diff(y(t),t,t)=yp4,diff(x(t),t)=Y[2],diff(y(t),t)=Y[4],x(t)=Y[1],y(t)=Y[3],{odes});
#The result of this is pasted as the equations 'eqs' in the procedure.
#YP[2] and YP[4] are found by fsolve.
p:=proc(N,t,Y,YP) local eqs,yp2,yp4;
   YP[1]:=Y[2];
   YP[3]:=Y[4];
   eqs:={yp2*sin(Y[4])+yp4^2 = -2*Y[1]*Y[3]*exp(-Y[2])+Y[1]*Y[4]*yp2,
        Y[1]*yp4*yp2+cos(yp4) = 3*Y[3]*Y[2]*exp(-Y[1])};
   YP[2],YP[4]:=op(subs(fsolve(eqs,{yp2,yp4}),[yp2,yp4]));
end proc:
res:=dsolve(numeric,procedure=p,initial=Array([1,0,0,1]),number=4,procvars=[x(t),diff(x(t),t),y(t),diff(y(t),t)],start=0);
plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..5);

Inside procedures you don't have full evaluation unlike the situation for commands from the input line.

That could be the problem, but if we see your definition of `*` we may be able to settle that question.

Your system is linear, so we could try using the tools available for such systems.

eqs:={-4.376768352*10^5*Mea-2.542105242*10^6*Agr-1.800849667*10^6*For-60170.80926*Res+8.115163696*10^6=0,-2.542105242*10^6*Mea-1.476500135*10^7*Agr-1.045965655*10^7*For-3.494828086*10^5*Res+4.713432037*10^7=0,-1.800849667*10^6*Mea-1.045965655*10^7*Agr-7.409712508*10^6*For-2.475766891*10^5*Res+3.339036627*10^7=0,-8272.145098*Res-60170.80926*Mea-3.494828086*10^5*Agr-2.475766891*10^5*For+1.115654125*10^6=0};
A,b:=LinearAlgebra:-GenerateMatrix(eqs,[Agr,For,Mea,Res]);
LinearAlgebra:-LinearSolve(A,b);
LinearAlgebra:-GaussianElimination(A);
LinearAlgebra:-ConditionNumber(A);

You will notice that the condition number is very high. Thus you can expect numerical problems. These are seen e.g. if you assign Digits to lower values than the default (10).

The classical methods use a fixed stepsize. In the help page

?dsolve,classical

we read:
'stepsize'
Specifies the static size of each step. Note that classical does not use error correction estimates to adapt the step size. The default step size is h = min((T0-Tf)/3, 0.005) for forward integration and h = max((T0-Tf)/3, -0.005) for backward integration, respectively, where the problem is solved over the interval from T0 to Tf

#You could try inserting a procedure which remembers what it computes. Here q(t) returns 1 if t is of type numeric and t otherwise.

restart:
q:=proc(t) option remember; if not type(t,numeric) then return 'procname(_passed)' else 1 end if end proc;
#Digits:= 32;
w :=  500;
ode := diff(y(t), `$`(t, 1)) = 2*I*y(t)*q(t) + sin(w*t)*y(t)*y(t):
ics := y(0) = 1:

p2 := dsolve({ics, ode}, numeric, method = classical[rk4],output=listprocedure,known=q):
T,Y:=op(subs(p2,[t,y(t)]));
#Finding the actual step:
T(1);
R:=op(4,eval(q)):
Lix:=remove(type,map(op,[indices(R)]),name):
L:=sort(Lix);
S:=seq(L[i]-L[i-1],i=2..nops(L));
{%};

#The reason you get 0.0025 and not 0.005 is that rk4 also computes the value of f(t,y(t)) in the middle of successive steps.

Given initial conditions at zero we need to determine alpha and beta so the conditions at t=3 are satisfied.

Here is one way:
restart;
odes := diff(x(t), t) = 4*x(t)-alpha*x(t)*y(t), diff(y(t), t) = -2*y(t)+beta*x(t)*y(t);
init := dsolve({odes, x(0) = 2, y(0) = 1}, numeric, parameters = [alpha, beta], output = listprocedure);
X, Y := op(subs(init, [x(t), y(t)]));
#Check
init(parameters = [1, 1]);
plot([X, Y], 0 .. 3);
X3:=proc(a,b)
   init(parameters=[a,b]);
   eval(X)(3)-4
end proc;
Y3:=proc(a,b)
   init(parameters=[a,b]);
   eval(Y)(3)-2
end proc;
res:=fsolve({X3,Y3},[1,1]);
init(parameters=res);
X(3),Y(3);
plot([X,Y],0..3);

with(LinearAlgebra,MatrixExponential);
A:=Matrix([[1,2,0],[0,1,2],[0,0,1]]);
f:=t-><2*t*exp(t),0,0>;
x0:=<2,1,0>;
X:=MatrixExponential(A,t).x0+int~(MatrixExponential(A,t-s).f(s),s=0..t);
#Check of solution.
eval(X,t=0);
diff~(X,t) = A.X+f(t);