Consider a different case:

M:=[1001, '$1..1000'$5000, 1001, '$1..1000'$5000, 1001]:  
time(ListTools:-SearchAll(1001, M));          
#                                    1.544
time(select(m->M[m]=1001, [seq](1..nops(M))));
#                                    14.628

So, it looks like SearchAll is attempting to be output sensitive and optimal in the case when there are very few occurrences. 

On the other hand, that makes this a O(n^2) worst case algorithm. Consider:

M:=['1001'$(10000)]:
time(ListTools:-SearchAll(1001, M));
#                                     0.992
M:=['1001'$(20000)]:                
time(ListTools:-SearchAll(1001, M));
#                                     3.924
M:=['1001'$(40000)]:                
time(ListTools:-SearchAll(1001, M));
#                                    15.996
time(select(m->M[m]=1001, [seq](1..nops(M))));
                                     0.044

Perhaps there should be flag to avoid the output sensitive case when you suspect there will be many occurences.

John

 It is doing a completely different computation:

sqrt(5*L[47]-L[48]): lprint(%);
3*8930^(1/2)

evalf(sqrt(5*L[47]-L[48]));
                                  283.4960317

John

 It is doing a completely different computation:

sqrt(5*L[47]-L[48]): lprint(%);
3*8930^(1/2)

evalf(sqrt(5*L[47]-L[48]));
                                  283.4960317

John

That a a great example of how to do the numerical computation very fast, but I should point out that it is not really a fair comparison.
For the exact operation this is about the same speed:

L:=[seq](2*i+1, i=10000..100000):
L:=Array(L):
L1:=ArrayTools:-Alias(L,[1..90000]):
L2:=ArrayTools:-Alias(L,1,[1..90000]):
time(map[inline](sqrt,5*L1+L2));
                                    2.708

John

That a a great example of how to do the numerical computation very fast, but I should point out that it is not really a fair comparison.
For the exact operation this is about the same speed:

L:=[seq](2*i+1, i=10000..100000):
L:=Array(L):
L1:=ArrayTools:-Alias(L,[1..90000]):
L2:=ArrayTools:-Alias(L,1,[1..90000]):
time(map[inline](sqrt,5*L1+L2));
                                    2.708

John

> L:=[seq](2*i+1, i=10000..100000):
> time(zip((x,y)->sqrt(5*x+y),L[1..-2],L[2..-1]));
                                     2.632
> restart
> L:=[seq](2*i+1, i=10000..100000):
> time(seq( sqrt(5*L[i]- L[i+1]),i=1..nops(L)-1) );
                                     1.676

(timings from a different machine than my earlier post)

John

> L:=[seq](2*i+1, i=10000..100000):
> time(zip((x,y)->sqrt(5*x+y),L[1..-2],L[2..-1]));
                                     2.632
> restart
> L:=[seq](2*i+1, i=10000..100000):
> time(seq( sqrt(5*L[i]- L[i+1]),i=1..nops(L)-1) );
                                     1.676

(timings from a different machine than my earlier post)

John

This is an automatic simplification that you can't do much about. Is it any comfort that it prettyprints nicely everywhere except prettyprint=0?

> a/(b*c);
                                       a
                                      ---
                                      b c
> lprint(%); # prettyprint=0 output 
a/b/c

> dismantle(%); # underlying DAG a^1*b^(-1)*c^(-1)

PROD(7)
   NAME(4): a
   INTPOS(2): 1
   NAME(4): b
   INTNEG(2): -1
   NAME(4): c
   INTNEG(2): -1

John

Interestingly, neither Maple 8 nor Maple 9 can solve this either. 

John

Interestingly, neither Maple 8 nor Maple 9 can solve this either. 

John

 Since _F1 is a function, you can substitute (almost) any function for it an have a solution.  The identity function is the simplest choice:
 

# sol is the output of pdsolve
id:=x->x:
sol1 := eval(sol, _F1=id):
expr := rhs(sol1);
                                          1
                          expr := --------------------
                                          3
                                 (y - 2 x)  (y + 2 x)

simplify(diff(expr,x)*(x+ y)+diff(expr,y)*(4*x + y));
                                       0

You could use use 'exp', or 'sin' instead of id, and it will also give a solution.  Anything sufficiently differentiable will work.  Your choice depends on your application I suppose.

John

 Since _F1 is a function, you can substitute (almost) any function for it an have a solution.  The identity function is the simplest choice:
 

# sol is the output of pdsolve
id:=x->x:
sol1 := eval(sol, _F1=id):
expr := rhs(sol1);
                                          1
                          expr := --------------------
                                          3
                                 (y - 2 x)  (y + 2 x)

simplify(diff(expr,x)*(x+ y)+diff(expr,y)*(4*x + y));
                                       0

You could use use 'exp', or 'sin' instead of id, and it will also give a solution.  Anything sufficiently differentiable will work.  Your choice depends on your application I suppose.

John

Very nice.

John

*edit: or that ---^

True, a one-liner would do it:

subsindets(sys, 'indexed', x->cat(op(0,x),'_',op(1,x),'_',op(2,x)) );

John

In fact it is a very subtle problem.  Indexed variables a[i,j] can cause problems in some solve subroutines, so one of the first things solve does is rename these variables to new symbols like x_nn.  In the case of this system, the order in which the polynomial system solver handles the variables matters a lot, and when the variables get renamed, they get ordered differently.  It seems more or less chance that the original variable names get ordered in a way so the polynomial system solver returns quickly.

Here is another work around (use symbols instead of indexed names):

vars := indets(sys);
vars2 := map(x->evaln(cat(op(0,x),'_',op(1,x),'_',op(2,x))), vars);
sys := subs([seq(vars[i]=vars2[i], i=1..nops(vars))], sys):
time(solve(sys));
                                    275.825

John