@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);

@emmantop Expanding exp(beta/theta) in a taylor series raises the question:
With respect to what variable and about what point?

Also as I was saying above even replacing exp(beta/theta) by 1 doesn't help.

Even trying replacing exp(...) by 1 as in

dsolve(eval(sys, {epsilon = E[1],exp=1}),numeric, method = bvp[midrich],maxmesh=256);

doesn't work.

When solving for the fourth derivative of f and expanding there is a term having
f''(eta)*theta'(eta) / ( f(eta)*theta(eta)^2)
and the exp-factor and a factor beta.

I tried setting beta=0 instead of replacing exp directly by 1. This has the effect of removing that term besides having the desirable side effect of setting exp=1.

Then there was no problem:

P1 := {Ec = .5, N = 10, Re = 5,lambda = -1, n = 2, sigma = 10, k[1] = .9}; #No beta value
#ode1, ode2, and bcs as before.
sys1:=eval([ode1, ode2, bcs], P1);
res:=dsolve(eval(sys1, {epsilon = E[1],beta=0}),numeric, method = bvp[midrich],maxmesh=256);
plots:-odeplot(res,[[eta,f(eta)],[eta,theta(eta)]],0..10);

####Looking at the term left out by setting beta=0 (for what it is worth):
T:=diff(theta(eta),eta)/theta(eta)^2;
F:=diff(f(eta),eta,eta)/f(eta);
plots:-odeplot(res,[[eta,F*T]],0..10)
#Huge at the end of the interval. 
 

@patient I detest document mode and 2D math. I have to fight it wasting valuable time. Personally I never use it.

Here is the worksheet version in 1D-math.

I added a dsolve-version at the end using events to stop integration a small amount before x = 2*x0 at which the exact solution for U(x) will become zero. Since U(x) appears in the denominator for V''(x) =  ... this will be critical for a numerical solution. Furthermore plotting is done there from x=0.001 in order to avoid the problem at x=0.

MaplePrimes15-05-09odesys_parametersC.mw

I have no solution, but notice that if you remove the exp-terms then there is no problem.

As in

dsolve(eval(sys, {epsilon = E[1],exp=0}),numeric, method = bvp[midrich],maxmesh=256);