Maybe it will be helpful to take a more general approach.  Your system has the form
diff(U,t) = diff(U,x$2) + A.U

where U = U(x,t) is a vector and A is a constant matrix.  If A can be diagonalized as
A = S.B.S^(-1) where B is diagonal, then the transformed vector V = S^(-1).U
satisfies the decoupled system

diff(V,t) = diff(V,x$2) + B.V

i.e. diff(v_j, t) = diff(v_j, x$2) + b_j * v_j

where b_j are the diagonal elements of B, the eigenvalues of A.  These in turn can be transformed into heat equations by

w_j(x,t) = exp(-b_j*t) * v_j(x,t)

These are "analytical" solutions.  There is no approximation involved.

Here's a version that returns a longest increasing sublist instead of just the length.

longestIncreasingSequence:= proc(L::list)
   local n,G,A,a,s,t,P,m,E,AE,q;
   uses GraphTheory; 
     n := nops(L);
     G := Digraph([$1..n,s,t],
                   { seq([[s,i],-1],i=1..n),
                      seq([[i,t],0],i=1..n),
                      seq(seq(`if`(L[i]<L[j]
                                  , [[i,j],-1]
                                  , NULL
                                 )
                             , j = i+1..n)
                         , i = 1..n-1)}
                   );
     A:= AllPairsDistance(G); 
     a:= s; 
     for m from A[n+1,n+2]+1 to 0 do
        E:= subs(t=NULL, Departures(G,a));
        AE:= map(v->A[v,n+2],E);
        member(m,AE,'q'); 
        a:= E[q]; P[m]:=a;
     end do;
     [seq](L[P[m]],m=A[n+1,n+2]+1 .. 0);
 end proc;

Here's a version that returns a longest increasing sublist instead of just the length.

longestIncreasingSequence:= proc(L::list)
   local n,G,A,a,s,t,P,m,E,AE,q;
   uses GraphTheory; 
     n := nops(L);
     G := Digraph([$1..n,s,t],
                   { seq([[s,i],-1],i=1..n),
                      seq([[i,t],0],i=1..n),
                      seq(seq(`if`(L[i]<L[j]
                                  , [[i,j],-1]
                                  , NULL
                                 )
                             , j = i+1..n)
                         , i = 1..n-1)}
                   );
     A:= AllPairsDistance(G); 
     a:= s; 
     for m from A[n+1,n+2]+1 to 0 do
        E:= subs(t=NULL, Departures(G,a));
        AE:= map(v->A[v,n+2],E);
        member(m,AE,'q'); 
        a:= E[q]; P[m]:=a;
     end do;
     [seq](L[P[m]],m=A[n+1,n+2]+1 .. 0);
 end proc;

pdsolve can't seem to handle this system, but I've found some solutions.

Look for solutions of the following form:

 Ansatz:= {u1(x,t)=A*exp(k1*x+k2*t), u2(x,t)=B*exp(k1*x+k2*t), 
                     u3(x,t)=C*exp(k1*x+k2*t)};
 des := {diff(u1(x,t),t)-(diff(diff(u1(x,t),x),x)
          -a*(u1(x,t)-u2(x,t))+b*(u3(x,t)-u1(x,t))),  
     diff(u2(x,t),t)-(diff(diff(u2(x,t),x),x)
          -c*(u2(x,t)-u3(x,t))+a*(u1(x,t)-u2(x,t))),
   diff(u3(x,t),t)-(diff(diff(u3(x,t),x),x)
          -b*(u3(x,t)-u1(x,t))+c*(u2(x,t)-u3(x,t)))};

See what equations have to be true for this to give a solution:

 eval(des, Ansatz);
 ABCeqs:= eval(%, exp=1); 
 with(LinearAlgebra):
 M := GenerateMatrix(ABCeqs, [A,B,C])[1];

                       [                          2]
                       [-b ,  -c , b + c + k2 - k1 ]
                       [                           ]
                  M := [                     2     ]
                       [-a ,  k2 + c + a - k1  , -c]
                       [                           ]
                       [   2                       ]
                       [-k1  + k2 + b + a , -a , -b]

For there to be a solution with A, B, C not all 0, this matrix must be singular.

 factor(Determinant(M));
 k2s := solve(%, {k2});

k2s := {k2 = k1^2}, 
  {k2 = -b-a-c+k1^2+(b^2-c*b-a*b+c^2-a*c+a^2)^(1/2)}, 
  {k2 = -b-a-c+k1^2-(b^2-c*b-a*b+c^2-a*c+a^2)^(1/2)}

for i from 1 to 3 do  
  N := op(NullSpace(eval(M, k2s[i])));
  sol[i] := eval(Ansatz, k2s[i] union {A = N[1], B = N[2],       
      C = N[3]})  
  end do; 

 You can then build solutions as linear combinations of sol[1], sol[2] and sol[3] for different k1.

Finding a maximum clique in a graph is NP-complete.  MaximumClique uses a branch-and-bound backtracking algorithm.  It should be simpler in this case because the digraph is transitive, but I don't know if MaximumClique would take advantage of this.  I tried a random list of length 100, and Joe's procedure took over 29 seconds.

Another way of looking at this is that you want the longest path in the digraph.  This would be the shortest path if you made each edge have "distance" -1.  The AllPairsDistance procedure finds the shortest distances between all pairs of vertices using the Floyd-Warshall algorithm (you can't use  Dijkstra's when there are negative distances).
 

longestIncreasingSequenceLength:= proc(L::list)
  local n,G,i,j,A;
  uses GraphTheory;
    n := nops(L);
    G := Digraph(n
                 , {seq(seq(`if`(L[i]<L[j]
                                 , [[i,j],-1]
                                 , NULL
                                )
                            , j = i+1..n)
                        , i = 1..n-1)}
                );
    1-min(op(map(rhs,rtable_elems(AllPairsDistance(G)))));
end proc;

This took just over 2 seconds on the same list of length 100.

Finding a maximum clique in a graph is NP-complete.  MaximumClique uses a branch-and-bound backtracking algorithm.  It should be simpler in this case because the digraph is transitive, but I don't know if MaximumClique would take advantage of this.  I tried a random list of length 100, and Joe's procedure took over 29 seconds.

Another way of looking at this is that you want the longest path in the digraph.  This would be the shortest path if you made each edge have "distance" -1.  The AllPairsDistance procedure finds the shortest distances between all pairs of vertices using the Floyd-Warshall algorithm (you can't use  Dijkstra's when there are negative distances).
 

longestIncreasingSequenceLength:= proc(L::list)
  local n,G,i,j,A;
  uses GraphTheory;
    n := nops(L);
    G := Digraph(n
                 , {seq(seq(`if`(L[i]<L[j]
                                 , [[i,j],-1]
                                 , NULL
                                )
                            , j = i+1..n)
                        , i = 1..n-1)}
                );
    1-min(op(map(rhs,rtable_elems(AllPairsDistance(G)))));
end proc;

This took just over 2 seconds on the same list of length 100.

It should return count+1 instead of count.

It should return count+1 instead of count.

Well, you could do this.  Note that you want the intersection of the regions

z - sqrt(x^2 + y^2) >= 0

z - 2 >= 0

4 - z >= 0

So this could be defined as min(z - sqrt(x^2+y^2), z - 2, 4 - z) >= 0.

Therefore you could say something like:

with(plots):
implicitplot3d(min(z - sqrt(x^2 + y^2), z - 2, 4 - z), 
  x = -5 .. 5, y = -5 .. 5, z = -5 .. 5);

 Or a little better:

P1:= implicitplot3d(min(z - sqrt(x^2 + y^2), z - 2, 4 - z),
    x = -4 .. 4, y = -4 .. 4, z = 1.5 .. 4.5,
   grid=[30,30,30],style=patchnogrid):
P2:= intersectplot(z=sqrt(x^2+y^2), z=2, x=-4 .. 4, y= -4 .. 4, z = 1.5 .. 2.5,
    colour=black):
P3:= intersectplot(z=sqrt(x^2+y^2), z=4, x=-4 .. 4, y = -4 .. 4, z = 3.5 .. 4.5,
   colour=black):
display(P1,P2,P3, axes=box, scaling=constrained);

Well, you could do this.  Note that you want the intersection of the regions

z - sqrt(x^2 + y^2) >= 0

z - 2 >= 0

4 - z >= 0

So this could be defined as min(z - sqrt(x^2+y^2), z - 2, 4 - z) >= 0.

Therefore you could say something like:

with(plots):
implicitplot3d(min(z - sqrt(x^2 + y^2), z - 2, 4 - z), 
  x = -5 .. 5, y = -5 .. 5, z = -5 .. 5);

 Or a little better:

P1:= implicitplot3d(min(z - sqrt(x^2 + y^2), z - 2, 4 - z),
    x = -4 .. 4, y = -4 .. 4, z = 1.5 .. 4.5,
   grid=[30,30,30],style=patchnogrid):
P2:= intersectplot(z=sqrt(x^2+y^2), z=2, x=-4 .. 4, y= -4 .. 4, z = 1.5 .. 2.5,
    colour=black):
P3:= intersectplot(z=sqrt(x^2+y^2), z=4, x=-4 .. 4, y = -4 .. 4, z = 3.5 .. 4.5,
   colour=black):
display(P1,P2,P3, axes=box, scaling=constrained);

It seems the help page ?solve,details is out of date and needs to be revised.  The Explicit option seems to have been introduced in Maple 10 without being documented.  Looking at the code for solve, other possible options are AllSolutions and split (I have no idea what that does).

Also, the statement about what _EnvExplicit does is incorrect.  Even with _EnvExplicit unassigned, Maple will give solutions in radicals for polynomials such as

x^9+3*x^7-3*x^6+3*x^5-6*x^4+5*x^3-3*x^2+4*x-1

whose roots lie in a field obtained from the rationals by cubic extensions
(in fact this polynomial is f(f(x)) + 2 where f(x) =  x^3 + x - 1).

Actually that's not quite right.  Everything entered at top level is evaluated, and that evaluation can cause commands to be executed. 

In the case of the document at hand, I don't think that's what's causing the slowness.  It seems to be the document mode and 2D math input.  On my (old and slow) computer, re-executing the worksheet typically took over 4 CPU seconds (as measured by time()).  A worksheet with the same commands in 1D input mode took only a small fraction of a second.

View 4541_da_cancellare.mw on MapleNet or Download 4541_da_cancellare.mw
View file details

Actually that's not quite right.  Everything entered at top level is evaluated, and that evaluation can cause commands to be executed. 

In the case of the document at hand, I don't think that's what's causing the slowness.  It seems to be the document mode and 2D math input.  On my (old and slow) computer, re-executing the worksheet typically took over 4 CPU seconds (as measured by time()).  A worksheet with the same commands in 1D input mode took only a small fraction of a second.

View 4541_da_cancellare.mw on MapleNet or Download 4541_da_cancellare.mw
View file details

Why did you say that wasn't what you wanted?  That was it exactly.  Or you could use implicitdiff to get an equivalent result.  Then replace diff(y(x),x) with -1/diff(y(x),x) to get a differential equation for the orthogonal trajectories.