@MotoVector Yes. After the command  g[5]:=89; you may try eval(g);

@J4James Maybe making a table with indices C,t and with a matrix containing the corresponding eta,u values  as entries would be useful:

restart;
Eq1:=diff(u(eta,t),t)=C*diff(u(eta,t),eta$2);
BCs := {u(0,t)=sin(t), u(L,t)=0};
ICs := {u(eta,0)=0};
L:=600:
sol:= C1->pdsolve(eval(Eq1,C=C1),ICs union BCs,numeric,time=t,range=0..L,spacestep=.1);
P:=(C1,t1)->proc(eta1) subs(sol(C1):-value()(eta1,t1),[eta,u(eta,t)]) end proc;
#I'm only making a smaller sized version, see the values of N, Cfinal, and tfinal
A:=table();
h:=.01: N:=9:
A[.01,.1]:=Matrix([seq(P(.01,.1)(h*i),i=0..N)]);
Cfinal:=.05:
tfinal:=0.4:
for C1 from 0.01 by 0.01 to Cfinal do
  for t1 from 0.1 by 0.1 to tfinal do
    A[C1,t1]:=Matrix([seq(P(.01,.1)(h*i),i=0..N)])
  end do
end do:

eval(A);
There may be faster ways of doing this.
A variant of P:
P:=(C1,t1)->proc(eta1) subs(sol(C1):-value(t=t1)(eta1),[eta,u(eta,t)]) end proc;


 

@J4James Maybe making a table with indices C,t and with a matrix containing the corresponding eta,u values  as entries would be useful:

restart;
Eq1:=diff(u(eta,t),t)=C*diff(u(eta,t),eta$2);
BCs := {u(0,t)=sin(t), u(L,t)=0};
ICs := {u(eta,0)=0};
L:=600:
sol:= C1->pdsolve(eval(Eq1,C=C1),ICs union BCs,numeric,time=t,range=0..L,spacestep=.1);
P:=(C1,t1)->proc(eta1) subs(sol(C1):-value()(eta1,t1),[eta,u(eta,t)]) end proc;
#I'm only making a smaller sized version, see the values of N, Cfinal, and tfinal
A:=table();
h:=.01: N:=9:
A[.01,.1]:=Matrix([seq(P(.01,.1)(h*i),i=0..N)]);
Cfinal:=.05:
tfinal:=0.4:
for C1 from 0.01 by 0.01 to Cfinal do
  for t1 from 0.1 by 0.1 to tfinal do
    A[C1,t1]:=Matrix([seq(P(.01,.1)(h*i),i=0..N)])
  end do
end do:

eval(A);
There may be faster ways of doing this.
A variant of P:
P:=(C1,t1)->proc(eta1) subs(sol(C1):-value(t=t1)(eta1),[eta,u(eta,t)]) end proc;


 

@Carl Love Maybe I'm not quite getting your point. You write:

"Suppose that I try to write a function with the same bahvaiour as `if`."

But then B and C ought to be evaluated explicitly inside the procedure as done below (?).

restart:
if2:= proc(A, B::uneval, C::uneval)
local r:= evalb(A);
     if r::truefalse then
          if r then eval(B) else eval(C) end if
     else
          'procname'(args)
     end if
end proc:
eval(if2(A, sin(Pi), C), A= true);
eval(`if`(A,sin(Pi),C), A=true);

@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