@Markiyan Hirnyk 

Sure:

CS:= L-> {L[]}:
L:= ['rand()' $ 2^20]:
CodeTools:-Usage(CS(L)):
memory used=8.00MiB, alloc change=0 bytes, cpu time=343.00ms, real time=357.00ms

Compare that with the near-neglible time for a full linear scan of L:

CodeTools:-Usage(member(-1,L)):
memory used=0.51KiB, alloc change=0 bytes, cpu time=16.00ms, real time=12.00ms

Are you saying that the text file literally contains the -6 on the line above the rest of the number? Or does the text file contain literally the characters "10^-6"?

@Markiyan Hirnyk 

This can also be done as

{L1[]} subset {L2[]};

However, such simplicity may obscure the fact that converting a long list to a set is a significant computation timewise.

@Joe Riel I hope the proprosal does not maintain the current limit of bases from 2 to 36. With the input in list form, the base can be anything.

Define what it means for one list to "contain" another. Does order matter? Does the repetition count of individual elements matter?

@Joe Riel I think that there is a library way to do it, however you may consider it cheating. See my answer below.

What does d* and drho mean? I presume some derivative? Where are your other ODEs? You'll need three ODEs.

You could just as well apply indices to a nested list:

PowerIthfactor:= (n::posint, i::posint)-> ifactors(n)[2,i,2];

It looks cleaner than op, IMO.

@ctc Let me know if you need a hint on the polynomial approximation.

Ke:= Matrix(9,9, (i,j)-> int(eval(Grad[i].Grad[j], xi= 0), eta= 0..1));

@diedrotn 

restart:

interface(showassumed= 0):
assume(xi::real, eta::real):
phi:= <
     (2*xi-1)*(xi-1)*(2*eta-1)*(eta-1),
     -4*xi*(xi-1)*(2*eta-1)*(eta-1),
     xi*(2*xi-1)*(2*eta-1)*(eta-1),
     -(8*xi-4)*(xi-1)*eta*(eta-1),
     16*xi*(xi-1)*eta*(eta-1),
     -4*xi*(2*xi-1)*eta*(eta-1),
     (2*xi-1)*(xi-1)*eta*(2*eta-1),
     -4*xi*(xi-1)*eta*(2*eta-1),
     xi*(2*xi-1)*eta*(2*eta-1)
>:

Gradient:= v-> <diff(v,xi), diff(v,eta)>:
Grad:= Gradient~(phi):
Ke:= Matrix(op(1,Grad), (i,j)-> int(Grad[i].Grad[j], [xi, eta]=~ 0..1));

It looks neater in the worksheet than it does here.

@Christopher2222 ln(1/x) = -ln(x) is only a rule for x > 0. Maple is solving over a larger domain. You can do this:

-ln(-a/(b-c));

combine(%, symbolic);

The uniqueness of the factorization can only be guaranteed when your factor base is pairwise relatively prime.

@nm Go to the Tools menu. Select Options. Select the first tab, General. Go to the second drop-down, "How should Maple handle the creation of new Math Engines?" Is your selection "Share one engine among all documents"? If yes, then change it to one of the other two options.

@Andriy 

When doing a type check with the function name in prefix form (e.g., as opposed to using specfunc), the number of anythings matches the number of arguments of the function.

type(f(x), f(anything, anything));
     false
type(f(x,y), f(anything, anything));     true
type(f(x,y,z), f(anything,anything));
     false