This variant covers the case where the first set has more than one element and also saves the result instead of just printing it.
It works for e.g.
L:=[{a}, {b, c}, {a, d, f}, {b, d}];
L:=[{a}, {b, c}, { d, f}, {b, d}];
L:=[{a,b}, {b, c}, { a,d, f,b}, {b, d}];

res:='res':
look:=L[1]:
for i from 2 to nops(L) do
 if look subset L[i] then res:=L[i] minus look; break end if
end do;
res;

@OffshoreEngineer Take a look at the following approach, where initially only the conditions at x = 0 are used.

W:=simplify(rhs(dsolve({deqns,inits[1..2]}))) assuming EL>0,IN>0;
indets(W,name) minus indets(deqns,name);
#We now need to impose the boundary conditions at x = 100.
IntegrationTools:-Expand(W);
W1:=collect(%,[_C1,_C3]); #Assuming that the arbitrary constants are indeed labelled 1 and 3
eq1:=eval(W1,x=100)=0;
eq2:=eval(diff(W1,x,x),x=100)=0;
Csol:=solve({eq1,eq2},{_C1,_C3}):
#The final solution is:
wsol:=eval(W1,Csol):
#However, notice that there are lots of integrals
indets(wsol,specfunc(anything,Int));
#Although they are all inert (i.e. Int instead of int) value doesn't help, since Maple cannot evaluate them, which this shows:
wsol-subs(int=Int,value(wsol)); #Result 0

Working with wsol will involve lots of numerical integration and it does not seem to be a good idea at all.

@OffshoreEngineer Take a look at the following approach, where initially only the conditions at x = 0 are used.

W:=simplify(rhs(dsolve({deqns,inits[1..2]}))) assuming EL>0,IN>0;
indets(W,name) minus indets(deqns,name);
#We now need to impose the boundary conditions at x = 100.
IntegrationTools:-Expand(W);
W1:=collect(%,[_C1,_C3]); #Assuming that the arbitrary constants are indeed labelled 1 and 3
eq1:=eval(W1,x=100)=0;
eq2:=eval(diff(W1,x,x),x=100)=0;
Csol:=solve({eq1,eq2},{_C1,_C3}):
#The final solution is:
wsol:=eval(W1,Csol):
#However, notice that there are lots of integrals
indets(wsol,specfunc(anything,Int));
#Although they are all inert (i.e. Int instead of int) value doesn't help, since Maple cannot evaluate them, which this shows:
wsol-subs(int=Int,value(wsol)); #Result 0

Working with wsol will involve lots of numerical integration and it does not seem to be a good idea at all.

@zippon I worked a little more on this.
I seem to be getting correct results all the time if I declare r local:
task := proc (a, i, j) local k,r;
etc.

@zippon I worked a little more on this.
I seem to be getting correct results all the time if I declare r local:
task := proc (a, i, j) local k,r;
etc.

Maybe you could provide examples. It is clearly not always a problem. Example:
A:=[seq](Int(x*exp(-x^(2*n)),x=-infinity..infinity),n=1..10);
value~(A);
evalf(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;

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