At closer inspection I now see that the angular frequencies appearing in sin and cos are v and 2*v, where
v = 0.5995250000.

I start as before with getting rid of that extra beta in beta(.787*beta). It should not have been there in the first place. I wonder how it got in there?

My guess is that the floats appearing hide a simple symbolic relationship beyond what I have already mentioned.
So we shall try to find that.

restart;
Digits:=10:
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = beta*(.787*beta) , A);
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);
indets(f,specfunc(anything,{sin,cos}));
subs(1.199050000*beta=2*v*beta,.5995250000*beta=v*beta,f);
op([2,1],%);
Cs:=coeffs(%,[beta,cos(2*beta*v),sin(2*beta*v),sin(beta*v)],'T');
T;
w:=add(a[i]*T[i],i=1..nops([T]));
param:= convert(zip(`=`,[seq(a[i],i=1..nops([T]))],[Cs]),set);
#For the integral to be convergent on the interval 0..1 it is necessary and sufficient that the series expansion:
series(w,beta=0,5);
#contains no powers of beta under 4:
eval(%,param union {v=.5995250000});

It so appears. But rounding as seen in my earlier response can very easily spoil the symbolic relationship:

a[2]+2*a[3]*v+a[4]*v^2 = 0

My advice is to keep symbolic input until you are forced to resort to numerics.

At closer inspection I now see that the angular frequencies appearing in sin and cos are v and 2*v, where
v = 0.5995250000.

I start as before with getting rid of that extra beta in beta(.787*beta). It should not have been there in the first place. I wonder how it got in there?

My guess is that the floats appearing hide a simple symbolic relationship beyond what I have already mentioned.
So we shall try to find that.

restart;
Digits:=10:
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = beta*(.787*beta) , A);
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);
indets(f,specfunc(anything,{sin,cos}));
subs(1.199050000*beta=2*v*beta,.5995250000*beta=v*beta,f);
op([2,1],%);
Cs:=coeffs(%,[beta,cos(2*beta*v),sin(2*beta*v),sin(beta*v)],'T');
T;
w:=add(a[i]*T[i],i=1..nops([T]));
param:= convert(zip(`=`,[seq(a[i],i=1..nops([T]))],[Cs]),set);
#For the integral to be convergent on the interval 0..1 it is necessary and sufficient that the series expansion:
series(w,beta=0,5);
#contains no powers of beta under 4:
eval(%,param union {v=.5995250000});

It so appears. But rounding as seen in my earlier response can very easily spoil the symbolic relationship:

a[2]+2*a[3]*v+a[4]*v^2 = 0

My advice is to keep symbolic input until you are forced to resort to numerics.

@mtj4u 

This time the function beta(.787*beta) appears. Does it mean (1) .787*beta, (2) beta, or (3) beta*(.787*beta)?
In any case f will have a singularity at zero sufficient to make the integral from 0 to 1 divergent.

restart;
Digits:=100:
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = beta*(.787*beta) , A); #First interpretation
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);
series(f, beta = 0, 9);

@mtj4u 

This time the function beta(.787*beta) appears. Does it mean (1) .787*beta, (2) beta, or (3) beta*(.787*beta)?
In any case f will have a singularity at zero sufficient to make the integral from 0 to 1 divergent.

restart;
Digits:=100:
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = beta*(.787*beta) , A); #First interpretation
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);
series(f, beta = 0, 9);

This time you somehow introduced another beta function: beta(.787*beta).
So the question is: what was the intention? Is it supposed to mean just  .787*beta, or is it supposed to mean beta, or beta*(.787*beta)?

In all 3 cases f has a serious singularity at zero making the integral over beta = 0.. 1 divergent.

restart;
Digits:=100: #overdoing it some maybe.
A := Int((1.334389725*sin(1.199050000*beta)/beta+(1.669308692*(cos(1.199050000*beta)-1.667987156*sin(1.199050000*beta)/beta+2.782181154*sin(.5995250000*beta)^2/beta^2))/beta^2)^2/((61.09*(240.38915*tanh(.787*beta)*beta(.787*beta)+1))/(61.09+3.935*tanh(.787*beta)*beta(.787*beta))+1.7314*coth(.787*beta)*beta(.787*beta)), beta = 0 .. infinity);

indets(A,specfunc(anything,beta));
subs(beta(.787*beta) = .787*beta, A); #Taking the first interpretation mentioned
B := normal(%);
f := IntegrationTools:-GetIntegrand(B);

asympt(f, beta, 9); #May not work in Maple 14 (but irrelevant)
series(f, beta = 0, 9);

My preference would be

evalindets(sol,specfunc(anything,{sin,cos}),s->op(0,s)(omega*t));

But your way could be made similarly:

evalindets(sol,specfunc(anything,{sin,cos}),x->applyrule(eq,x));

My preference would be

evalindets(sol,specfunc(anything,{sin,cos}),s->op(0,s)(omega*t));

But your way could be made similarly:

evalindets(sol,specfunc(anything,{sin,cos}),x->applyrule(eq,x));

@mapleq2013 x(t)-1 is a rootfinding trigger, meaning that when during integration x(t) - 1 hits zero the action (in this case 'halt') is activated.
So in the example there are two events both of the type simple rootfinding and both having action 'halt'.

Here is a simple example of a conditional trigger:

sys2:=diff(x(t),t)= y(t) , diff(y(t),t)= -x(t);  
res2:=dsolve({sys2,x(0)=1/2,y(0)=0},numeric,events=[[[x(t),y(t)>0],halt]]);
plots:-odeplot(res2,[t,x(t)],0..5);

In this one integration stops when x(t) = 0 provided that also y(t) > 0.

@mapleq2013 x(t)-1 is a rootfinding trigger, meaning that when during integration x(t) - 1 hits zero the action (in this case 'halt') is activated.
So in the example there are two events both of the type simple rootfinding and both having action 'halt'.

Here is a simple example of a conditional trigger:

sys2:=diff(x(t),t)= y(t) , diff(y(t),t)= -x(t);  
res2:=dsolve({sys2,x(0)=1/2,y(0)=0},numeric,events=[[[x(t),y(t)>0],halt]]);
plots:-odeplot(res2,[t,x(t)],0..5);

In this one integration stops when x(t) = 0 provided that also y(t) > 0.

@mtj4u The result from the integration over the finite interval r..infinity is the result of the integral from 0 to infinity of f with an error less than the value of the sum of epsilon and the integral of g from r to infinity. In the code epsilon = 1e-10 and really means relative error, but since the result is close to 1, there is no difference between relative and absolute error. 

Since I observed some error messages in your worksheet I upload my version with r = 10^6.

MaplePrimesIntegral1.mw

Edited: 'sum of' has replaced 'maximum of'. The text in the worksheet has been adjusted accordingly.

NB. You are using Maple 14 it appears. I'm using Maple 16. The maxintervals option for the methods _d01ajc and _d01akc is not available in Maple 12 I just found out. I have no access to Maple 14, but that option is available in Maple 15.

@mtj4u The result from the integration over the finite interval r..infinity is the result of the integral from 0 to infinity of f with an error less than the value of the sum of epsilon and the integral of g from r to infinity. In the code epsilon = 1e-10 and really means relative error, but since the result is close to 1, there is no difference between relative and absolute error. 

Since I observed some error messages in your worksheet I upload my version with r = 10^6.

MaplePrimesIntegral1.mw

Edited: 'sum of' has replaced 'maximum of'. The text in the worksheet has been adjusted accordingly.

NB. You are using Maple 14 it appears. I'm using Maple 16. The maxintervals option for the methods _d01ajc and _d01akc is not available in Maple 12 I just found out. I have no access to Maple 14, but that option is available in Maple 15.

@PlpPlp I changed things a little to reflect my understanding of the problem. I use worksheet interface and 1D input (so much easier for me).

DeltaF := K[1]*m[1]^(2)*m[2]^(2)+(K[1]+(B[1]^(2))/(2* c[11])-(B[2]^(2))/(2 *c[44]))*(m[1]^(2)+m[2]^(2))*m[3]^(2)+K[2]*m[1]^(2)*m[2]^(2)*m[3]^(2)+B[1]*(u[m1]*m[1]^(2)+u[m2]*m[2]^(2))+B[2]*u[m6]*m[1]^(2)*m[2]^(2)-B[1]*((B[1])/(6 *c[11])+(c[12])/(c[11])*(u[m1]+u[m2]))*m[3]^(2)+1/(2)*mu[0]*M[s]^(2)*(N[1]*m[1]^(2)+N[2]*m[2]^(2)+N[3]*m[3]^(2))-mu[0]*M[s]*(H[1]*m[1]+H[2]*m[2]+H[3]*m[3]);

H[eff1] := -diff(DeltaF, m[1])/mu[0];
H[eff2] := -diff(DeltaF, m[2])/mu[0];
H[eff3] := -diff(DeltaF, m[3])/mu[0];

expr := collect(mu[0]*(m[2]*H[eff3]-m[3]*H[eff2]),[m[1],m[2],m[3]],distributed,factor);

#Now linearising about m = m0. The number 2 at the end means that the orders neglected are 2 and above.

mtaylor(expr, [m[1]=m[1,0],m[2]=m[2,0],m[3]=m[3,0] ], 2);

@PlpPlp I changed things a little to reflect my understanding of the problem. I use worksheet interface and 1D input (so much easier for me).

DeltaF := K[1]*m[1]^(2)*m[2]^(2)+(K[1]+(B[1]^(2))/(2* c[11])-(B[2]^(2))/(2 *c[44]))*(m[1]^(2)+m[2]^(2))*m[3]^(2)+K[2]*m[1]^(2)*m[2]^(2)*m[3]^(2)+B[1]*(u[m1]*m[1]^(2)+u[m2]*m[2]^(2))+B[2]*u[m6]*m[1]^(2)*m[2]^(2)-B[1]*((B[1])/(6 *c[11])+(c[12])/(c[11])*(u[m1]+u[m2]))*m[3]^(2)+1/(2)*mu[0]*M[s]^(2)*(N[1]*m[1]^(2)+N[2]*m[2]^(2)+N[3]*m[3]^(2))-mu[0]*M[s]*(H[1]*m[1]+H[2]*m[2]+H[3]*m[3]);

H[eff1] := -diff(DeltaF, m[1])/mu[0];
H[eff2] := -diff(DeltaF, m[2])/mu[0];
H[eff3] := -diff(DeltaF, m[3])/mu[0];

expr := collect(mu[0]*(m[2]*H[eff3]-m[3]*H[eff2]),[m[1],m[2],m[3]],distributed,factor);

#Now linearising about m = m0. The number 2 at the end means that the orders neglected are 2 and above.

mtaylor(expr, [m[1]=m[1,0],m[2]=m[2,0],m[3]=m[3,0] ], 2);

You have posted several of those by now. However, I fail to see the connection to mathematical software like Maple.

Maybe uploading the worksheet with your example could help make the problem clear.