@J4James Here is an ordered version (revised from an earlier version to make it faster):

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);
hC:=.01: ht:=.1: he:=.01: #Increments in C,t, and eta, respectively
Ne:=4: NC:=3: Nt:=2: #Number of increments
A:=Matrix(NC*Nt*Ne,4); #Rows: [C,eta,t,u]
t0:=time():
r:=0:
for i from 1 to NC do
   q:=sol(i*hC);
   for j from 1 to Nt do
     v:=q:-value(t=j*ht);
     for k from 1 to Ne do
       r:=r+1;
       A[r,..]:=subs(v(he*k),Vector[row]([i*hC,eta,t,u(eta,t)]))
     end do
   end do
end do;
time()-t0;
#interface(rtablesize=infinity);
A;

@J4James Here is an ordered version (revised from an earlier version to make it faster):

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);
hC:=.01: ht:=.1: he:=.01: #Increments in C,t, and eta, respectively
Ne:=4: NC:=3: Nt:=2: #Number of increments
A:=Matrix(NC*Nt*Ne,4); #Rows: [C,eta,t,u]
t0:=time():
r:=0:
for i from 1 to NC do
   q:=sol(i*hC);
   for j from 1 to Nt do
     v:=q:-value(t=j*ht);
     for k from 1 to Ne do
       r:=r+1;
       A[r,..]:=subs(v(he*k),Vector[row]([i*hC,eta,t,u(eta,t)]))
     end do
   end do
end do;
time()-t0;
#interface(rtablesize=infinity);
A;

@geri23 It may be worth your while to formulate your problem again from the very start and then post that as a new question in MaplePrimes. That way other people get a chance to participate in the discussion. The present thread is likely read by only you and me.

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