Unless you explicity assigned mf as an Array and are using parentheses as indices, the assignment you are doing is probably not what you want. Regardless, it isn't clear what you want.  If you just want to create a list of items corresponding to the rhs of your equation, then do

L := [seq(seq(f(Prc)*m(Prc,B), B = 2..200,0.5), Prc = 2..100)];
sort(L);

Generally, for digit comparison, it is useful to convert to strings and compare characters.  Here the ?StringTools[CommonPrefix] procedure should be useful.

I believe my hypothesis can be fixed by constraining the applicable functions.  Here's the basic idea

restart;
with(IntegrationTools):

F := 1/deltax*Int(sin(f(x)),x=x1..x1+deltax);
Ftheta1 := Change(F, f(x)=theta, [theta]);
Ftheta2 := eval(Ftheta1, f(x1+deltax)=convert(series(f(x1+deltax),deltax,2),'polynom'));
Ftheta3 := eval(Ftheta2, [f(x1)=2*n*Pi, deltax = Pi/D(f)(x1)]);
Fphi1 := Change(Ftheta3, phi=theta-2*n*Pi, [phi]);
Fphi2 := eval(Fphi1, sin(phi+2*n*Pi) = sin(phi));
Fphi3 := evalindets(Fphi2
                    , specfunc(specfunc(anything,RootOf),D(f))
                    , DF -> eval(series( DF, phi, 2), RootOf(-f(_Z) + 2*n*Pi) = x1)
                   );

                   );
                           Pi
                          /
               D(f)(x1)  |                   sin(phi)
      Fphi3 := --------  |    -------------------------------------- dphi
                  Pi     |                 (2)
                         |               (D   )(f)(x1)            2
                        /     D(f)(x1) + ------------- phi + O(phi )
                          0                D(f)(x1)

This suggests that if f''(x)/f'(x)^2 approaches 0 as x goes to infinity, then the limit goes to 2/Pi.  That explains why ln(x) fails.  However, both x^2 and x*exp(x) meet the critieria.  One can directly compute the limit for f = ln(x):

restart;
with(IntegrationTools):
f := ln:
F := 1/deltax*Int(sin(f(x)),x=x1..x1+deltax);
Ftheta := Change(F, f(x)=theta, [theta]);
sol1 := solve(f(x1)=2*n*Pi,[x1]);
sol2 := solve(f(x1+deltax)=(2*n+1)*Pi,[deltax]);
y1 := simplify(eval(eval(Ftheta,sol2[1]),sol1[1]),ln) assuming real;
y2 := Change(y1, theta0=theta-2*n*Pi, [theta0]);
y3 := eval(y2, sin(theta0+2*n*Pi)=sin(theta0));
y4 := value(y3);
y5 := simplify(expand(y4));
                                        1 + exp(Pi)
                             y5 := 1/2 -----------
                                       exp(Pi) - 1
                

Seems like you might be able to handle the digits by using logarithms.  Euclid would then be changed to compute the sum of the logs.  Regardless, this problem quickly becomes intractable.

To determine where the slow parts are, create a procedure and then profile it.  You can use the ?CodeTools[Profiling] package. 

restart;
euclid := proc (s, n)
    if 1 < s and 1 <= n then
        mul(1/(1-ithprime(i)^(-s)), i = 1 .. n);  # changed product to mul; it is faster
    end if
end proc:

Print := proc()
local q,m,s,n,b;
    for q to 5 do       # reduced to 5 to speed up the run
        m := 10^(-q);
        for s from 2 to 10 do
            n := 1;
            b := 1;
            while b >= m do
                b := evalf(abs(Zeta(s)-euclid(s, n)));
                n := n+1
            end do;
            #print(m, s, n-1)    # remove output for now
        end do
    end do
end proc:

with(CodeTools:-Profiling):

t := Build('procs'=[euclid, ithprime, Zeta, Print], 'commands'='Print()'):
PrintProfiles(t);
SortBy(t, order=time);

function                 calls      time     time%          words    words%
---------------------------------------------------------------------------
euclid                    2288     8.547     57.90      166746250     85.84
ithprime               1686699     5.386     36.49       13812972      7.11
Print                        1     0.812      5.50       13643772      7.02
Zeta                       129     0.016      0.11          40111      0.02
---------------------------------------------------------------------------
total:                 1689117    14.761    100.00      194243105    100.00

Clearly euclid is the problem.  There is a simple way to significantly improve its performance as used; create a recursive  version that caches its results.  That works because euclid is repeatedly called with increasing values of n, so we can reuse the previous result.

euclid := proc(s,n)
option cache;
    if n = 1 then
        1/(1-1/2^s);
    else
        procname(s,n-1)/(1-1/ithprime(n)^s);
    end if;
end proc:

Reprofiling gives

SortBy(t, order=time);
function                 calls      time     time%          words    words%
---------------------------------------------------------------------------
Print                        1     0.828     88.09       13628634     79.07
ithprime                  5941     0.056      5.96         634754      3.68
euclid                    1963     0.044      4.68        2935657     17.03
Zeta                       129     0.012      1.28          37970      0.22
---------------------------------------------------------------------------
total:                    8034     0.940    100.00       17237015    100.00

I expect both to be 2/Pi. 

F(x) := whatever;

is not equivalent to

F := x -> whatever;

Both are legal, but the first is rarely what you want. It doesn't create a procedure, rather it merely assigns whatever to the function call F(x), which has almost nothing to do with, say, F(3).  As an example, consider

F := x -> x:    # assign a procedure
F(3) := 100: # assign an entry in F's remember table

We can now evaluate

F(z);
            z
F(3);
            100

What has actually happened is that the procedure F is assigned but also has a remember table that bypasses the computation when the argument is precisely 3.  Note that the argument must match exactly for the entry in the remember table to be accessed. For example,

F(3.) := 24:
F(3), F(3.), F(3.0);
              100, 24, 3.0

For your example, what you probably want is

F := unapply(int(x^2, x = 0..x), x);

Clearer might be to use separate variables for the independent variable of the integral and its upper limit:

F := unapply(int(x^2, x=0..y), y);

I was asking what dsys was assigned, you left that out of the steps shown. Probably it is

dsys := { ec1, ec2, <initial conditions> };

where <initial conditions> are just that, maybe alpha(0) = 0, beta(0) = 0. 

A minor problem is that tau[1] and tau[2] are not assigned numeric values.  They could be declared as parameters to dsolve, but at some point they would have to get numeric values before integrating. The more serious issue is that dsolve cannot convert your set of differential equations to a first order system.

The usual way to combine plots is with ?plots[display].  For example,

y1 := x+4:
y2 := 2*x+1:
p1 := plot(y1, x=0..5, color=red):
p2 := plot(y2, x=0..10, color=blue):
plots:-display(p1,p2);

Don't worry about the code appearing to be "cut off"; we can copy the entire block without a problem.

A nice way to solve the problem is to extend your system so that the integration is performed by dsolve.  To do that, just add the differential equation diff(tot(t),t) = S(t) + R(t), with tot(0) = 0.  Thus,

with(plots):

deqs := {diff(S(t), t) = t^2
        , diff(R(t), t) = t
        , diff(tot(t),t) = S(t)+R(t)
       }:

ini := {S(0) = 1, R(0) = 1, tot(0)=0}:

ode1 := dsolve(deqs union ini, numeric
               , 'events' = [[t = 600, S(t) = 10^7], [t = 800, S(t) = 5*10^7]]):

odeplot(ode1, [t, S(t)+R(t), color = blue, thickness = 2], 0 .. 1000, numpoints = 1000); 
ode1(1000);
[t = 1000., R(t) = 500001.000000000, S(t) = 0.212666666666667 10 ,

                                 11
    tot(t) = 0.469666682666667 10  ]

What do you plan to do with them?  You might consider writing a simple module that provides some conversions.  For example,

krylov := module()
option package;

global K1,K2,K3,K4;

export To, From, Simplify, Evalf;

local ModuleLoad, removeImaginary;

    ModuleLoad := proc()
    global `convert/tokrylov`, `convert/fromkrylov`, `simplify/krylov`
        , `evalf/K1`, `evalf/K2`, `evalf/K3`, `evalf/K4`;

        `convert/tokrylov`   := To;
        `convert/fromkrylov` := From;
        `simplify/krylov`    := Simplify;
        `evalf/K1`           := x -> Evalf(K1(x));
        `evalf/K2`           := x -> Evalf(K2(x));
        `evalf/K3`           := x -> Evalf(K3(x));
        `evalf/K4`           := x -> Evalf(K3(x));

    end proc;

    To := proc(ex)
        eval(ex, [cos = K1-K3,
                  sin = K2-K4,
                  cosh = K1+K3,
                  sinh = K2+K4
                 ]);
    end proc:

    From := proc(ex)
        eval(ex, [K1 = 1/2*cosh + 1/2*cos,
                  K2 = 1/2*sinh + 1/2*sin,
                  K3 = 1/2*cosh - 1/2*cos,
                  K4 = 1/2*sinh - 1/2*sin
                 ]);
    end proc:

    Evalf := proc(ex)
        evalf(From(ex), _rest);
    end proc;

    Simplify := proc(ex)
        # a very basic simplification, Kx(I*z) --> Kx(z)
        evalindets(ex
                   , 'specfunc(anything,{K1,K2,K3,K4})'
                   , fx -> op(0,fx)(removeImaginary(op(1,ex)))
                  );
    end proc:

    removeImaginary := proc(ex)
        if ex :: imaginary
        or ex :: `*` and ormap(type,ex,imaginary)
        then
            return ex/I;
        else
            return ex;
        end if;
    end proc:

    ModuleLoad();

end module:

with(krylov):

Simplify(K1(I*x));
                                     K1(x)

convert(exp(2*x),'compose,trigh,tokrylov');
                     K1(2 x) + K3(2 x) + K2(2 x) + K4(2 x)

evalf(K1(Pi/2 + K3(3)),30);
                        303.220352664355942804059979113

This is just something I threw together in a few minutes; you'd want to rewrite it to do something useful

You need to reformulate the equation for a(t) so that it is not in terms of a(t).  The dsolve command would normally do that, but the piecewise hinders it. You also need to replace time with t. The following should do what you want

m := 5; F0 := 10; a0 := 5;
deqs := { m*a(t) = piecewise(F0/m*t < a0, F0*t, 0)
          , v(t) = diff(x(t), t)
          , a(t) = diff(v(t), t)
        }:
ics := {x(0) = 0, v(0) = 0};
integ := dsolve( deqs union ics, numeric):

plots:-odeplot(integ, [t,x(t)], 0..3);

Note that for the units to be consistent, F0 has dimensions of force/time.

Followup

I see that does not work past t=2.5.  One way around that is to use dsolve event handling. While this can be a bit complicated, in this case we can handle it simply by using a discrete variable (delta) whose value changes from 1 to 0 when F0/m*t-a crosses zero. A better model might pay attention to the direction---not an issue here with the given equations and initial conditions.

deqs := { m*a(t) = F0*t*delta(t)
          , v(t) = diff(x(t), t)
          , a(t) = diff(v(t), t)
        }:

ics := {x(0) = 0, v(0) = 0, delta(0)=1};

# Use a single event that is triggered when the acceleration
# reaches the threshold.  At that time, change delta to 0
# which causes the effective acceleration to go to 0.
evts := [[F0/m*t-a0=0, delta(t)=0]]:

integ := dsolve( deqs union ics, 'numeric'
                 , 'events' = evts
                 , 'discrete_variables' = [delta(t)::float]
               ):

plots:-odeplot(integ, [t,v(t)], 0..3);

The error message indicates what needs to be done.  In the assignment to P, change each x to x(t).  Alas, dsolve won't give a symbolic solution, but you can use dsolve,numeric.

integ := dsolve({eq,ics}, numeric):
integ(1);
plots:-odeplot(integ, [t,x(t)], 0..1);

@Chef-Chaudard It would be helpful if you showed both what you did and what Maple returned.  Are you sure you replaced the single forward quotes ('...') with double quotes ("...")? The double quote is a single symbol, not two single quotes.

The format operand to fprintf should be a string; a name may work but only if Maple can parse it, and here it won't parse the $, which is an operator. So enclose the format operand in double quotes.