Remember also to replace x by X e.g. like this:
P:=subs(x=X,g1)+g22-R*(g33+g44):
and also in eta:
Eta:=subs(x=X,eta);
P1:=(1/Eta)*Int(P,y=0..Eta):

Remember also to replace x by X e.g. like this:
P:=subs(x=X,g1)+g22-R*(g33+g44):
and also in eta:
Eta:=subs(x=X,eta);
P1:=(1/Eta)*Int(P,y=0..Eta):

@emma hassan 

A := Matrix([seq(`ϕ`[1, i],i=1..m)]);

@emma hassan 

A := Matrix([seq(`ϕ`[1, i],i=1..m)]);

@Carl Love No, you are right. fsolve doesn't seem to care at all. Initially I introduced the weights for plotting purposes and with the intent of minimizing the weighted sum of squares. With eps = 1e-3 the contributions were of comparable size.
(The procedure name p2 was ill chosen as that is also the name of the parameter in eq2.)

@Carl Love No, you are right. fsolve doesn't seem to care at all. Initially I introduced the weights for plotting purposes and with the intent of minimizing the weighted sum of squares. With eps = 1e-3 the contributions were of comparable size.
(The procedure name p2 was ill chosen as that is also the name of the parameter in eq2.)

@samiyare The following seems to work OK.
I'm applying different weights to the two outputs (eps and eps^(-1)).

p2:=proc(pp2,phi0,eps::numeric:=1e-3) option remember; local res,F0,F1,F2,c1,c2;
#print([_passed]); #Remove # to see progress
res := dsolve({eq1=0,subs(p2=pp2,eq2)=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=phi0,T(0)=0,D(T)(0)=1}, numeric,output=listprocedure);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
c1:=evalf(2/(1-zet^2)*Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*( F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )/(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet))-pp2;
c2:=evalf(2/(1-zet^2)*(Int((1-eta)*F1(eta),eta=0..1-zet)))-0.02;
c1*eps,c2/eps
end proc;
p2a:=proc() p2(_passed)[1] end proc;
p2b:=proc() p2(_passed)[2] end proc;
S:=(x,y)->ln(p2a(x,y)^2+p2b(x,y)^2): #For plotting
plot3d(S,61000..63000,0.007..0.009,grid=[15,15],view=-5..4);
res:=fsolve([p2a,p2b],[61800..62200,0.0082..0.0084]);
#Returns [61966.095945736166, 0.82962774008845973e-2]
#And
p2(op(res));
#returns -1.2*10^(-16), 7.000000000000000000*10^(-17)

@samiyare The following seems to work OK.
I'm applying different weights to the two outputs (eps and eps^(-1)).

p2:=proc(pp2,phi0,eps::numeric:=1e-3) option remember; local res,F0,F1,F2,c1,c2;
#print([_passed]); #Remove # to see progress
res := dsolve({eq1=0,subs(p2=pp2,eq2)=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=phi0,T(0)=0,D(T)(0)=1}, numeric,output=listprocedure);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
c1:=evalf(2/(1-zet^2)*Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*( F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )/(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet))-pp2;
c2:=evalf(2/(1-zet^2)*(Int((1-eta)*F1(eta),eta=0..1-zet)))-0.02;
c1*eps,c2/eps
end proc;
p2a:=proc() p2(_passed)[1] end proc;
p2b:=proc() p2(_passed)[2] end proc;
S:=(x,y)->ln(p2a(x,y)^2+p2b(x,y)^2): #For plotting
plot3d(S,61000..63000,0.007..0.009,grid=[15,15],view=-5..4);
res:=fsolve([p2a,p2b],[61800..62200,0.0082..0.0084]);
#Returns [61966.095945736166, 0.82962774008845973e-2]
#And
p2(op(res));
#returns -1.2*10^(-16), 7.000000000000000000*10^(-17)

@Carl Love After a little plotting:

fsolve(p(pp2)=0,pp2=63000..63500);
                          63181.95053

@Carl Love After a little plotting:

fsolve(p(pp2)=0,pp2=63000..63500);
                          63181.95053

This comment belongs further down, but for some reason I'm not able to put it there on the little machine I'm using.
You could set the 306 free b's equal to zero (or any other number).
Using the notation in your worksheet:
seq(evalb(VVV(t)[i,1]-AAA[i,1]=0),i=1..18);
(VVV(t)[3,1]-AAA[3,1]);
indets(K1(t));
op(0,b[1,2]);
bb:=select(x->evalb(op(0,x)=b),indets(K1(t)));
bb0:=zip(`=`,convert(bb,list),[0$nops(bb)]);
K2:=subs(bb0,K1(t));
indets(K2);

Since the solution remains the same (without giving the method) if the initial conditions are all zero, i.e.

ICS := x1(0) = 0, x2(0) = 0, D(x1)(0) = 0, D(x2)(0) = 0;
it seems quite fair to call this a bug.
But at least a warning should be issued maybe saying something like "Gave up trying to satisfy the initial conditions".

Since the solution remains the same (without giving the method) if the initial conditions are all zero, i.e.

ICS := x1(0) = 0, x2(0) = 0, D(x1)(0) = 0, D(x2)(0) = 0;
it seems quite fair to call this a bug.
But at least a warning should be issued maybe saying something like "Gave up trying to satisfy the initial conditions".

@awass As for your third point it seems that the parameters version is faster and less costly.
You could try (in both cases) the following:
interface(warnlevel=0):
CodeTools:-Usage(seq(P(k*.01),k=1..10000)):
For the parameters version I got
memory used=378.97MiB, alloc change=4.00MiB, cpu time=14.38s, real time=15.37s
For the version in your point 3 I got
memory used=1.33GiB, alloc change=36.00MiB, cpu time=42.34s, real time=43.24s

@awass As for your third point it seems that the parameters version is faster and less costly.
You could try (in both cases) the following:
interface(warnlevel=0):
CodeTools:-Usage(seq(P(k*.01),k=1..10000)):
For the parameters version I got
memory used=378.97MiB, alloc change=4.00MiB, cpu time=14.38s, real time=15.37s
For the version in your point 3 I got
memory used=1.33GiB, alloc change=36.00MiB, cpu time=42.34s, real time=43.24s