Since q, w, and e are procedures, I don't understand what you mean by using them as the variables with respect to which the integrations are done.

Do you not understand the error message? The differential order is 6: 3 from the Eq3 and 2 for theta and 1 for f from the Eq4. You have seven conditions. Eliminate 1 of them.

@Markiyan Hirnyk 

It's not big, but it is superlinear. But I'm not saying that it is a bad thing to do! I'm just saying that one has to watch out for seemingly trivial things which, if done in an inner loop, can add up. Note that converting a set to a list is relatively fast (but I can't figure out why it is not O(1).)

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