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

@patient I'm afraid we have returned to one of the problems I pointed out in my original answer: You cannot give initial conditions at z=0, which here corresponds to x=0.

Another unrelated point:
To check the exact solution given somewhat implicitly by SOL you need tro isolate U(x) and V(x):

solve(SOL,{U(x),V(x)});
odetest(%,SYS);
SOL2:=simplify(%) assuming real;

##############################
I now choose initial conditions taken from the exact solution.
I use z=-1 (arbitrarily: We need z<>0).
eval(SOL2,{x=x0,t=-x0^2});
ICS1:=simplify(%) assuming x0>0;
convert(diff(SOL2[2],x),D);
eval(%,x=x0);
ICS2:=eval(%,{t=-x0^2});
ICS:=ICS1 union {ICS2};
SYS2:=simplify(subs(t=-x0^2,SYS)) assuming x0>0;
numer(lhs(SYS2[2]))/x0^5;

res:=dsolve(SYS2 union ICS,numeric,parameters=[x0],method=rkf45_dae,maxfun=100000);
p := proc (t, itvl::range) res(parameters = [sqrt(-t)]); plots:-odeplot(res, [[x, U(x)], [x, V(x)]], itvl, _rest) end proc;
plots:-animate(p,[t,0..2,thickness=3],t=-1..-.1); p1:=%:
plots:-animate(plot,[subs(SOL2,[U(x),V(x)]),x=0..2,color=[red,green]],t=-1..-0.1); p2:=%:
plots:-display(p1,p2);

@javaneh0 After making the corrections that are pointed out by Makiyan then using method=laplace will give you a shorter looking result:

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

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

vars:=indets({sys_ode},function(identical(t)));
dsolve({sys_ode,ics},vars,method=laplace);

To get a shorter or more concrete result you need to give concrete values to the parameters, i.e. to the names yc, yh, etc.
You may in that case consider solving numerically.
Example of solving numerically:
All parameters are set to the same value 1 (easy to do that is why):

params:=indets({sys_ode},name) minus {t}=~1;
dsolve({sys_ode,ics},vars,method=laplace);
res:=dsolve(eval({sys_ode,ics},params),numeric);
plots:-odeplot(res,[t,ph(t)],0..5);