What you have in the first integral is essentially this:

Int(x^j*exp(-a*x),x=0..infinity);
value(%) assuming a>0,j>-1;
simplify(%) assuming j::posint; #If that is indeed the case?
convert(%,factorial);

Try this:

eq1,eq2:=diff(f1(x),x)=0.25*b*(a/b-f1(x))^2, diff(f2(x),x)=0.25*b*(a/b-f2(x)+f1(x))^2;
ic1,ic2:=f1(0)=0,f2(0)=0;
dsolve({eq1,eq2,ic1,ic2});
simplify(%);
#The (escaped local) names _Z1~ and   _Z3~ (or similar) represent arbitrary integers:
indets(%,`local`);
If you put them equal to zero you get the same result as when you solve in two steps as done afterwards:
res1:=eval(%%,%=~0);
dsolve({eq1,ic1});
subs(%,{eq2,ic2});
res2:=dsolve(%);
subs(res1,f2(x))-subs(res2,f2(x));
simplify(%);
#In Maple 16 it simplifies to 0.

#Comment: If you do instead res1:=eval(%%,%=~1); you get a result for which f2(0) <> 0. Since the initial value has a unique solution we know in advance that either the solution doesn't in fact depend on the escaped locals (the _Z's) or the solution is wrong for some sets of _Z's.

You need to replace G by G(t,z[1],z[2],z[3]), but I don't see much of a chance for obtaining a solution with Maple:

Eq1 := subs(  G=G(t,z[1],z[2],z[3]),'((z[1]-1)*lambda+(z[3]-1)*gamma)*G+(1-z[1])*delta*(diff(G, z[1]))+(N*z[3]-z[2])*kappa*(diff(G, z[2]))+(1-z[3])*mu*(diff(G, z[3]))+(z[2]-z[1]*z[3])*beta*(diff(diff(G, z[3]), z[1])) = diff(G, t)' );
#pdsolve separates variables, but that is it:
pdsolve(Eq1);

Try this:

A:=Int(Int((1-x)*(1-y)*(x+y)*(x^n_1)*(y^n_2)*(1-x-y)^n_3, y=0..1-x),x=0..1);
#I have used Int instead of int so we see the integral before proceeding with a calculation.
#Decent looking results (here in Maple 12) are only available for concrete values of the exponents:
eval(A,{n_1=1,n_2=2,n_3=3});
value(%);
#But you may try
value(A) assuming posint;

Not many corrections were needed.
restart;
#A, B = Payoff matrices for players
A := Matrix(2, 2, [4, 2, 1, 3]); B := Matrix(2, 2, [1, 4, 3, 2]);
#T defines the time interval [0,T], s is the step length
T := 80; s := 1/5; k := [seq(i*s, i = 1 .. T)];
#Trajectories for different initial conditions
S1 := dsolve({z1(0) = .5, z2(0) = .18, diff(z1(t), t) = z1(t)*(1-z1(t))*(A[1, 2]-A[2, 2]-(A[1, 2]+A[2, 1]-A[1, 1]-A[2, 2])*z2(t)), diff(z2(t), t) = z2(t)*(1-z2(t))*(B[1, 2]-B[2, 2]-(B[1, 2]+B[2, 1]-B[1, 1]-B[2, 2])*z1(t))}, numeric,output=Array(k));
S2 := dsolve({z1(0) = .5, z2(0) = .1, diff(z1(t), t) = z1(t)*(1-z1(t))*(A[1, 2]-A[2, 2]-(A[1, 2]+A[2, 1]-A[1, 1]-A[2, 2])*z2(t)), diff(z2(t), t) = z2(t)*(1-z2(t))*(B[1, 2]-B[2, 2]-(B[1, 2]+B[2, 1]-B[1, 1]-B[2, 2])*z1(t))}, numeric,output=Array(k));
#Trajectories
g1 := [seq([S1[2, 1][i, 3], S1[2, 1][i, 2]], i = 1 .. T)]:
g2 := [seq([S2[2, 1][i, 3], S2[2, 1][i, 2]], i = 1 .. T)]:
plot([g1, g2], axes = boxed, color = [black,blue]);

I don't see any immediate solution, but the following at least works.

C:=Matrix([[D[1],0,D[2]],[0,D[2],D[1]]]);
B:=;
freeze~(B);
C.%;
evalindets(%,`*`,apply@op);
thaw~(%);
%;

#Putting this into a procedure:

`&p`:=proc(d,v) local a,b;
   a:=d.freeze~(v);
   b:=thaw~(evalindets(a,`*`,apply@op));
   eval(b)
end proc;
   
C &p B;

#####Later addition:
In the more general case where your differential operators in the matrix are polynomials in D[1] and D[2] you will have to freeze C too:

Take as a simple example the matrix A below which has the linear differential operator T as one of its elements:

T:=-3*D[1]+7*D[2]+(f->6*f);
#Testing it on a function g:
T(g);
A:=Matrix([[T,0,D[2]],[0,D[2],D[1]]]);
#Now just modify the procedure above so that A will be frozen too:
`&p`:=proc(d,v) local a,b;
   a:=freeze~(d).freeze~(v);
   b:=thaw~(evalindets(a,`*`,apply@op));
   eval(b)
end proc;
   
A &p B;

I have corrected typos and whatever to get the following version, which works.

restart:

Eq1:=diff(f(eta),eta,eta,eta)+0.5*f(eta)*diff(f(eta),eta,eta)=0;

Eq2:=diff(theta(eta),eta,eta)+0.5*Pr*f(eta)*diff(theta(eta),eta)=0;

bc1:=f(0)=0, D(f)(0)=0, D(f)(10)=1;

#No []
fixedparameter:=Pr=0.72;

#eval not evalf
Eq3:=eval(Eq2,fixedparameter);

#a*... ? and theta(eta) probably should be theta(0) or ?
bc2:=theta(10)=0, D(theta)(0)= -a*(1-theta(0));

L:=[0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1, 5, 10, 20];

for k from 1 to 10 do

#bc1,bc2 (no s)
R:=dsolve(eval({Eq1, Eq3, bc1, bc2},

a= L[k]), [f(eta), theta(eta)], numeric, output = listprocedure);


Y||k :=rhs(R[5]); Yp||k := -rhs(R[6]);

end do:

#plot, not plots
plot([Y||(1..10)],0..6,labels= [eta, thet(eta)]);

[Yp||(1...10)(0)];
[Y||(1..10)(0)];

(Comment: I changed the names of the procedure for my own use. That of course is irrelevant.)
Using the result from LSSolve with Le = 8 as an initialpoint for LSSolve when Le=10 I was able to get a result:

With Le=5 first do

res:=Optimization:-LSSolve([pPhi,pTheta,pF1,pS]);
inipt:=convert(res[2],list);
#Keep inipt for later use
#Next set Le = 8 and do
Optimization:-LSSolve([pPhi,pTheta,pF1,pS],initialpoint=inipt);
#Again keep the result for later use in inipt2, say.
#Finally set Le = 10 and do
Optimization:-LSSolve([pPhi,pTheta,pF1,pS],initialpoint=inipt2);

#You can get a result for Le = 8 without option 'initialpoint', but that requires adding the option maxfun=0 in the call to dsolve/numeric. It takes a while then.

#Added July 5:
Here is a faster version which uses option remember in a procedure 'paramset', whose only job is to set the parameters.
Since only the procedural version of LSSolve is used there is no need for the lines about returning unevaluated on nonnumeric input.

restart;
eq1:=diff(f(x), x$3)+diff(f(x), x$2)*(A*f(x)+x/2)+diff(f(x), x$1)*(1- A*diff(f(x),x$1))+A-1  ;
eq2:=diff(s(x), x$2)+diff(s(x), x$1)*(A*f(x)+x/2)+s(x)*(1-2*diff(f(x),x$1));
eq3:=diff(theta(x), x$2)/Pr+A*( f(x)*diff(theta(x),x$1))+x*0.5*diff(theta(x), x$1)+Nb*diff(theta(x),x$1)*diff(phi(x),x$1)+Nt*diff(theta(x),x$1)^2;
eq4:=diff(phi(x), x$2)+(A*f(x)+x*0.5)*diff(phi(x), x$1)*Le+Nt*diff(theta(x),x$2)/Nb;

b := 8:  A:=2: Nb:=0.1: Nt:=0.1: Pr:=1:
#I let 'Le' be unassigned, but include it as a parameter in dsolve.
#Le:=8:
solpar:=dsolve({eq1=0,eq2=0,eq3=0,eq4=0, f(0)=0, D(f)(0)=0, D(D(f))(0)=f2,s(0)=1, D(s)(0)=s1,theta(0)=1, D(theta)(0)=t1,phi(0)=1, D(phi)(0)=p1}, parameters=[Le,f2,t1,p1,s1],numeric,output=listprocedure,maxfun=0):
phin,thetan,sn,f1n:=op(subs(solpar,[phi(x),theta(x),s(x),diff(f(x),x)])):

paramset:=proc(Le,f2,t1,c1,s1) option remember; solpar(parameters=[Le,f2,t1,c1,s1]); NULL end proc;
#pPhi computes phi(b)
pPhi:=proc(f2,t1,p1,s1)
paramset(Le,_passed);
phin(b)
end proc:
#pTheta computes theta(b)
pTheta:=proc(f2,t1,p1,s1)
paramset(Le,_passed);
thetan(b)
end proc:

#pS computes sn(b)
pS:=proc(f2,t1,p1,s1)
paramset(Le,_passed);
sn(b)
end proc:

#pF1 computes f1n(b)-1
pF1:=proc(f2,t1,p1,s1)
paramset(Le,_passed);
f1n(b)-1
end proc:

forget(paramset); #Erases the memory table just in case you do this several times
t0:=time():
inipt:=table():
inipt[4]:=[]:
for Le from 5 to 10 do
  res:=Optimization:-LSSolve([pPhi,pTheta,pF1,pS],initialpoint=inipt[Le-1]);
  inipt[Le]:=convert(res[2],list)
end do;
time()-t0;

#There is an interesting problem when Le=9.

Try

indets({Solseq1},`local`);
indets({Solseq2},`local`);

The copied version will after execution refer to the global names d and `a[1]`.

restart;
f:=(-31250000 + 312500 *(x + y) +
 625 *(4* a* x - 4* x^2 - 5* x* y + 4 *(b - y)* y) +
 25* x* (-a* b + b* (x - y) + y^2))/((-500 + 5* a)* (500 - 5* b));
g:=(-1/(1200* (100 - a)* (-100 + b)))*(-2 *a^3* (-100 + b) + 200* b^3 - b^4 +
   6* a* (-100 + b)* (200 - x)* x + 120000* x^2 - 800* x^3 + x^4 +
   120000* y^2 - 1200* x* y^2 - 800* y^3 + 8 *x* y^3 +
   6* b* (x^3 - 20000* y - x* y^2 + 100* (-2* x^2 + 2* x* y + y^2)));
#I hope I got your f and g right.
pxy:=proc(aa,bb) global f,x,y; local F,res;
F:=subs(a=aa,b=bb,f);
res:=Optimization:-Minimize(F,x=aa..100,y=bb..100);
op(subs(res[2],[x,y]))
end proc;
#Testing:
pxy(.4,5);
G:=unapply(g,a,b,x,y);
gg:=(a,b)-> G(a,b,pxy(a,b));
Testing:
gg(1,2);
Optimization:-Maximize(gg,0..100,0..100);

Does the following simple example suit your situation. The original equation is simply diff(x(t),t)=(-x(t)^2-1).

sys:=diff(x(t),t)=(-x(t)^2-1)*piecewise(x(t)>=0,1,0);
sol:=dsolve({sys,x(0)=1},numeric);
plots:-odeplot(sol,[t,x(t)],0..10,thickness=2);

You can solve eq1 with the initial conditions first and then use the 'known' option:

sol1:=dsolve({eq1 =0,f(0) = S, D(f)(0) = 1, (D@D)(f)(0) = slop}, numeric,output=listprocedure);
F:=subs(sol1,f(x));

sol:=dsolve(subs(f(x)=F(x),{eq2=0,eq3=0, C(0)=1,C(b)=0, T(0)=1, T(b) = 0}), numeric,output=listprocedure,known=F);

sol(0);

odeplot(sol1,[x,diff(f(x),x)],0..b);
odeplot(sol,[x,T(x)],0..b);
odeplot(sol,[x,C(x)],0..b);

In fact there is no solution. Assume to the contrary that F is a solution. Observing (as in Hartman, Ordinary Differential Equations 1964, p. 532) that
F''' + 27*F*F'' = 0 is a first order linear differential equation in the unknown F'' we obtain

(*) F''(x) = C*exp(-27*int(F(t),t=0..x)).

Now by the requirement F'(-infinity) = 0 the constant cannot be negative since then F' would be decreasing. C cannot be zero either, since then F' would be constant. Both violate the other requirement F'(infinity) = 1. Thus C > 0, which means that F' is increasing, thus in particular F'(x) > 0 for all x. Since F(-infinity) = 0 it follows that F(x) > 0 for all x.

Integrating (*) from x1 to 0 we get

(**) F'(0) - F'(x1) = C*int( exp(-27*int(F(t),t=0..x)), x=x1..0).

The left hand side of (**) has a limit as x1 -> -infinity (by assumption).
However, exp(-27*int(F(t),t=0..x)) = exp(27*int(F(t),t=x..0)) > exp(0) = 1 for all x < 0.
Thus the right hand side of (**) doesn't have a (finite) limit as x1 -> -infinity. So we have a contradiction.

The two seemingly different results are in fact equal. Try 'combine' on the difference.

Notice that if the boundary value problem has a solution, then it has infinitely many.

Suppose F were a solution to the bvp. Then also G(x) = F(x+p) would be a solution to the bvp for any constant p.
Since F cannot be periodic because of the boundary values imposed, there would be infinitely many solutions.