@acer eval evaluates its argument fully like most other procedures. The behavior of eval when the fully evaluated first argument does not contain the lhs of the second argument (i.e. r) seems reasonable. Why evaluate when the argument is already fully evaluated?
Since `eval/if` is using eval, I suppose the behavior is not a bug.
restart;
f:=x->x;
eval('f'(2),r=9); #Full evaluation of 'f'(2) is f(2), which has no r.
eval('f'(2,r),r=9);#Full evaluation of 'f'(2,r) is f(2,r), which has an r.
showstat(`eval/if`);
`eval/if`(`if`(r,f(2),p),r=true); #`if`doesn't evaluate its 2. and 3. arguments before calling `evalf/if`
#Same output as from
`eval/if`([r,'f(2)',p],r=true);

@acer The appearance of r in the arguments to f has an effect on eval:
restart;
f:=x->x;
eval('f'(2),r=9);
eval('f'(2,r),r=9);
eval(`if`(r,f(2+r),p),r=true);
eval(`if`(r,f(2,r),p),r=true);

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".