@Kitonum Certainly this is much nicer. I don't know what the original question was hinting at, but if it had to do with exhibiting a weakness in verify then it succeeded.

@Harry Garst In the help page for verify/equal it says:

The verify(expr1, expr2, equal) function returns true if
                signum(0, expr1 - expr2, 0) = 0
 and false if the call to signum returns a numeric object. Otherwise, it returns FAIL.


If you try

signum(0,abs(x)- sqrt(x^2), 0);

in Maple 2015 you get 
i.e. you don't get 0 nor a numeric object. Thus verify must return FAIL according to its help page.
However, if you use simplify on the signum-result then you get 0.
###################
I haven't been using verify (ever?) and cannot say that I'm impressed:
verify(I*sin(7),cos(8),equal);
signum(0,I*sin(7)-cos(8), 0);
verify(I*sin(7.),cos(8.),equal);
signum(0,I*sin(7.)-cos(8.), 0);

@Athar Shahabinejad I think I understood what you meant. The problem is that you cannot find a formula for the eigenvalues and thus not for the eigenvectors either. Therefore you cannot find a formula for the solution to your odes. Hence no formula for pf(t) and consequently you cannot evaluate the integral.

@Christopher2222 You cannot have an event written like you do:

[diff(theta(t),t)=0,a:=[op(a),t]]

The assignment is clearly out of place, but a=[op(a),t] would also be wrong.
Before I try to elaborate I should like to know what 'a' is supposed to be. What has it got to do with EL. Is it supposed to be a discrete variable?

My guess is that you intended something like

sol3 := dsolve({EL, ini} union {b(0)=[]},numeric,discrete_variables=[b(t)::list],events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)],[diff(theta(t),t)=0,b(t)=[op(b(t)),theta(t)]]],range=0..10);

But you get an error telling you that discrete variables must be of type boolean, integer, or float.

This restriction is explicitly mentioned in the help page
?dsolve,numeric,events
under the heading Discrete variables:

Only the types boolean, integer and float are currently available.

################
If my guess is right you may try something like this instead:

sol3 := dsolve({EL, ini} union {b(0)=0,TT(0)=0},numeric,discrete_variables=[b(t)::float,TT(t)::float],events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)],[diff(theta(t),t)=0,[TT(t)=t,b(t)=theta(t)]]],output=listprocedure);
B,TTp:=op(subs(sol3,[b(t),TT(t)]));
plots:-odeplot(sol3,[[t,theta(t)],[t,b(t)],[t,TT(t)]],0..10);
L:=[]:
for tt from 0 to 10 by .5 do L:=[op(L),[TTp(tt),B(tt)]] end do:
L;
ListTools:-MakeUnique(L);

@Christopher2222 How about
solL := dsolve({EL, ini},numeric,events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)]],range=0..10,output=listprocedure);
TH,TH1:=op(subs(solL,[theta(t),diff(theta(t),t)]));
fsolve(TH1,[2]);
TH(%);
fsolve(TH1,[4]);
TH(%);
## or
sol2 := dsolve({EL, ini} union {b(0)=0},numeric,discrete_variables=[b(t)::float],events=[[theta(t)=0,diff(theta(t),t)=-diff(theta(t),t)],[[diff(theta(t),t)=0,theta(t)>0],b(t)=theta(t)]],range=0..10);
plots:-odeplot(sol2,[[t,theta(t)],[t,b(t)]],0..10);
sol2(5);

@Athar Shahabinejad I see no way of doing that.
The characteristic polynomial is a polynomial of degree 6 with rather general coefficients. In fact it appears that they are so general that any monic polynomial of degree 6 can be written in that form. Thus it is well known that no formula exists (or can be found) for the 6 roots of that polynomial.
Unless you give values to all or some of the parameters I cannot see that you can get anywhere.

@Athar Shahabinejad
When I download directly from MaplePrimes my own worksheet, which I gave a link to above, open it in Maple 2015 and execute the whole thing in its entirety, I get that the integral is 7200.67239251219.

What did you do?

@Athar Shahabinejad OK, now the equation

eq:=ph(t)+pc(t)+pa(t)+pgd(t)+prd(t)+prj(t)+pf(t)=1  for all t

is satisfied because of the odes and the initial conditions. Thus it is not an extra requirement.
However, since we are using floats and consequently approximations, there will be small errors in the computation of any of the 7 functions. This is of no real interest when plotting, but may affect the computation of the integral
int(1-pf(t), t=0..infinity).

In fact almost surely it will, since we shall find that pf(t) -> almost exactly 1, but not quite 1, as t-> infinity.
Any deviation from 1 is disastrous to convergence.

Thus I shall advocate using the second approach I gave above, which meant eliminating pf using the equation eq above.
MaplePrimes15-05-15odeLINEAR_F.mw

############# int versus numerical integration ##############
If you use numerical integration then it turns out that that is clever enough to handle the problem I described.
Suppose you use the first approach I gave (i.e. not eliminating pf). Then both of these work:
evalf(Int(add(SOL[i],i=1..6),t=0..infinity));
evalf(Int(1-SOL[7],t=0..infinity));
but neither of these do:
int(add(SOL[i],i=1..6),t=0..infinity);
int(1-SOL[7],t=0..infinity);

@Athar Shahabinejad With the revised sys_ode I get as output from
simplify(add(u,u=sys_ode));

that the derivative of the sum of dependent variables is 2*prj(t), i.e. not identically zero.

So although the system is different, the divergence of the integral in question still holds.

@Athar Shahabinejad Although I don't know the background for the problem I have the sneaking suspicion that your odes are set up incorrectly; that in fact the derivative of the sum of all the dependent variables should be identically zero.
So I suggest checking sys_ode. Ought it not be that
simplify(add(u,u=sys_ode));
returns a zero on the right hand side, meaning that the sum of the dependent variables is constant and therefore with initial conditions in mind equal to 1 for all t.
The situation reminds me of radioactive decay A -> B -> C -> ... -> Z (stable), where the total number of atoms is constant and where it will surely also follow from the odes.

@Athar Shahabinejad In that case you need to eliminate one of the odes since the 7 odes together with the initial conditions already determine the solution uniquely.
In other words: unless the requirement 
ph(t)+pc(t)+pa(t)+pgd(t)+prd(t)+prj(t)+pf(t)=1
is automatically satisfied adding this extra requirement will mean that there won't be any solution.
If you add all 7 odes you will notice that the derivative of the sum above is not identically zero (see below). Thus you must abandon one ode. I chose the last one.

restart;
Digits:=15:
sys_ode := diff(ph(t),t) = (1-yc)*pc(t)+yh*prj(t)+urd*prd(t)+ugd*pgd(t)-yc*ph(t), diff(pc(t),t) = yc*ph(t)-(2-yc)*pc(t), diff(pa(t),t) = ya*pc(t)-pf(t), diff(prj(t),t) = yrj*pa(t)+prj(t), diff(prd(t),t) = -urd*prd(t)+yrd*pa(t), diff(pgd(t),t) = -ugd*pgd(t)+ygd*pa(t), diff(pf(t),t) = (1-ygd-yrj-yrd)*pa(t)+(1-yh)*prj(t);
simplify(add(u,u=sys_ode)); #Right hand side not identically zero!!!
ics := ph(0) = 1, pc(0) = 0, pa(0) = 0, prj(0) = 0, prd(0) = 0, pgd(0) = 0, pf(0) = 0;
eq:=ph(t)+pc(t)+pa(t)+pgd(t)+prd(t)+prj(t)+pf(t)=1;
pfres:=solve(eq,{pf(t)});
sys2:=eval([sys_ode[1..6]],pfres);
VARS:=map2(op,1,lhs~(sys2));
A,b:=LinearAlgebra:-GenerateMatrix(rhs~(sys2),VARS);
#The system is now inhomogeneous:
diff~(Vector(VARS),t)=A.Vector(VARS)-b;
p:=LinearAlgebra:-CharacteristicPolynomial(A,lambda);
params:={yc=1/3600,yh=1,yrd=1/180,ygd=1/180,yrj=1/180,urd=0.8,ugd=0.8,ya=1/180};
eval(p,params);
fsolve(%=0,lambda,complex);
L,V:=LinearAlgebra:-Eigenvectors(evalf(eval(A,params)));
VARS;
ics2:=ics[1..6]; #Leaving out info about pf
ICS:=Vector(rhs~([ics2]));
#A particular solution:
solp:=evalf(eval(A,params))^(-1).b;
#The solution
sol:=solp+V.LinearAlgebra:-DiagonalMatrix(exp~(t*L)).V^(-1).(ICS-solp);
evalc(sol);
SOL:=simplify(fnormal~(%));
plot(SOL,t=0..15);
SOL;
odetest(VARS=~SOL,eval(sys2,params));
fnormal(%);
simplify(%);
pfres;
one_minus_pf:=add(u,u=SOL);
int(one_minus_pf,t=0..infinity);
#Again infinity

@
Just two observations:

I let your worksheet open in the standard GUI. There I had to remove the `&epsilon` (or whatever it is) and write plain epsilon in the dsolve commands (not in the odes).

Changing theta(N) = 0 to D(theta)(N)=0 avoids the exp-problem since theta(eta) doesn't get too near zero (at least for N=2). 

@javaneh0 You say that you forgot the initial condition "ph+pc+pa+pgd+prd+prj+pf=1".

Since you are setting
ics := ph(0) = 1, pc(0) = 0, pa(0) = 0, prj(0) = 0, prd(0) = 0, pgd(0) = 0, pf(0) = 0;

this condition is not additional, it is automatically satisfied.

Maybe you mean the requirement ph(t)+pc(t)+pa(t)+pgd(t)+prd(t)+prj(t)+pf(t)=1  for all t?

@javaneh0 My code above is plain text and can be copied and pasted into a fresh Maple worksheet without any problem.

If you have problems with that let me know.

@patient Sorry about venting my frustration with the document mode and 2D-math. It is not your fault. It is the default environment in Maple these days if you don't change it, which, however, is easily done under Tools/Options.

For the exact solution just use piecewise:
solp:={u(z)=piecewise(z>-4,(1+z/4)*exp(z/8),0),v(z)=piecewise(z>-4,exp(z/4), -4*exp(-1)*1/z)};
SOLp := subs(u(z) = (-t)^(3/2)*U(x), v(z) = -V(x)*t, z = x^2/t, solp);

SOL2p:=solve(SOLp,{U(x),V(x)});
odetest(SOL2p,SYS);
simplify(%) assuming x^2/t>-4;
simplify(%%) assuming x^2/t<-4;

plots:-animate(plot,[subs(SOL2p,[U(x),V(x)]),x=0..2,color=[red,green],thickness=3],t=-1..-0.1);
For the numerical solution I see no reason not to just use the exact solution after integration stops.
After all it is a somewhat artificial stop that is made. The exact solution for U is continuous, yes, but not differentiable at the point where it becomes zero.

resE:=dsolve(SYS2 union ICS,numeric,parameters=[x0],method=rkf45_dae,events=[[x=2*x0-.001,halt]]);
pE := proc (t, itvl::range) local p1,p2;
  resE(parameters = [sqrt(-t)]);
  p1:=plots:-odeplot(resE, [[x, U(x)], [x, V(x)]], itvl, _rest);
  p2:=plot([0,4*exp(-1)/x^2],x=2*sqrt(-t)..rhs(itvl),_rest);
  plots:-display(p1,p2)
end proc;
plots:-animate(pE,[t,0.001..2,thickness=3],t=-1..-.1);