Removing the absolute values (which are not necessary for real solutions), the system (without initial conditions) has solutions involving integrals of RootOf's of complicated algebraic integrals:

> sys:= {diff(x[1](t),t$2) = k[1]/(x[2](t)-x[1](t))^2,
                 diff(x[2](t),t$2) = k[2]/(x[2](t)-x[1](t))^2});
   dsolve(sys);

What this comes from, basically, is a differential equation obtained by eliminating the variable x[2].

> DEtools[casesplit](sys);

That fourth-order differential equation for x[1](t) depends only on the second and higher derivatives of x[1](t), so we can think of it as a second-order equation in a = (D@@2)(x[1]).

> eval(op([1,2],%),diff(x[1](t),t$2)=a(t));

And this has implicit solutions involving one integration of an algebraic function.

> dsolve(%);

Anyway, all this isn't going to help you much if you're looking for explicit solutions satisfying your initial conditions.  For particular values of the parameters k[1] and k[2], you can get numerical solutions.

> S:= dsolve({eval(sys_ode,{k[1]=1,k[2]=2}), ICS}, numeric);
    plots[odeplot](S,[[t,x[1](t)],[t,x[2](t)]], t= 0 .. 20);

Is this what you mean?

> R3values:= [0,-1,-3,-4]; colours:= [red,blue,green,black];
   display([seq(plottools[transform](
      unapply([k, subs(parvalues, Pec_i = -20, R3=R3values[i], F2)], (k,S)))
      (implicitplot(subs(parvalues,Theta0=4,R3=R3values[i],Pec_i=-20,A1),
              k=0..20, S=0..1, colour=colours[i],gridrefine=3)), i=1..4)], 
      labels=[k,'F2']);

>asympt(RootOf(8*x^3+y^3-6*x*y-3,y),x);

Note that for any x there will be three complex y values (counted by multiplicity).  For large x these will correspond to the three values of (-1)^(2/3).  Presumably what you want is the real one, so take (-1)^(1/3) = -1 and (-1)^(2/3) = 1.

With symbolic quantities Q and beta in your expressions, which ones have positive imaginary parts is likely to depend on the values of those quantities.  However, you
can try finding those where the imaginary part is positive at certain values of Q and beta. 


> select(t -> is(evalc(Im(subs(beta=1,Q=1,t)))>0), t1);

Expanding on hirnyk's suggestion, here's how (in Standard GUI) you can get %f(2) to look like 2^2.

>   f:= x -> x^2;
     `print/f`:= x -> Typesetting:-msup(Typesetting:-Typeset(x),Typesetting:-mn("2"));

Then

> %f(2);
   value(%);

You might try this one. 

> labelledpolygon:= proc(L::list,r::positive,col)
# L is a list of 2D points, possibly including labels as the third elements
# r the distance from the vertices to the labels
# col the colour for the polygon
  uses plots;
  local V, T,P1,V2,i,a,b,c,n;
  V:= map(t -> [t[1],t[2]],L);
  if nops(L[1]) = 3 then T:= map(t -> t[3],L)
  else T:= map(convert,[$1..nops(L)],string)
  end if; 
  P1:= pointplot([op(V),V[1]],style=line,colour=col);
  n:= nops(V); b:= V[1]-V[n]; b:= evalf(b/sqrt(b[1]^2+b[2]^2));
  for i from 1 to n do
    a:= b;
    if i = n then b:= V[1]-V[n] else b:= V[i+1]-V[i] end if;
    b:= evalf(b/sqrt(b[1]^2+b[2]^2));
    c:= evalc([Re,Im](-(a[1]+I*a[2])*sqrt(-(b[1]+I*b[2])/(a[1]+I*a[2]))));
    V2[i]:= [op(V[i]+r*c),T[i]];
  end do;
  display(P1,textplot(convert(V2,list)),scaling=constrained,_rest); 
end proc;

Then for example:

> labelledpolygon([[1,2,"a"],[2,5,"b"],[4,0,"c"]],0.5, blue);

> labelledpolygon([seq([cos(i*Pi/4),sin(i*Pi/4)],i=0..7)],0.2,red,axes=box);

You might use the parametric form of dsolve(..., numeric).  See the help page ?dsolve,numeric,interactive.  For putting several solution curves (produced by odeplot) in the same plot, use display from the plots package.

As Christopher222 mentioned, building a list one element at a time is a bad idea. 
Building a sequence one element at a time is not really better.  But a Vector can be built efficiently one element at a time, even if you don't know in advance how many elements it should have, using "programmer indexing" with parentheses () rather than brackets [].
See the help page ?rtable_indexing. 

Compare:

> restart; ti:= time():
   S:= NULL:
   for i from 1 to 10000 do S:= S,i od:
   L:= [S]:
   time()-ti;
   
        5.803

> restart;
   ti:= time():
   V:= Vector():
   for i from 1 to 10000 do V(i):= i od:
   time()-ti;   

          0.110

I do not know of any way to modify the palettes.

Your pde is

> diff(h(x, t), t) = -.5*(3*s*(diff(h(x, t), x, x, x, x))-3*s*(diff(h(x, t), x, x))/a^2-A*s*(12*(diff(h(x, t), x))^2/h(x, t)^5-3*(diff(h(x, t), x, x))/h(x, t)^4+1749.6*s^5*(diff(h(x, t), x, x))/(h(x, t)^10*d^6)-17496*s^5*(diff(h(x, t), x))^2/(h(x, t)^11*d^6))*(2*h(x, t)^3*(1/3))/(12*Pi*h0^3))+3*(3*s*(diff(h(x, t), x, x, x))-3*s*(diff(h(x, t), x))/a^2-A*s*(-3*(diff(h(x, t), x))/h(x, t)^4+1749.6*s^5*(diff(h(x, t), x))/(h(x, t)^10*d^6))/(12*Pi*h0^3))*(diff(h(x, t), x))*h(x, t)^2;

Note: it's much better to post text rather than images of equations.  Fortunately I could recover the "Associated Text" from the image info.

It's very unlikely that such a complicated nonlinear PDE will have useful closed-form solutions (except in this case you might note that constants are solutions).  So you'll probably want to solve it numerically.  In order to do that, you must specify numerical values for the parameters, as well as appropriate initial and boundary conditions.  Then pdsolve(..., numeric) should do it for you.

For example:

> pde1:= eval%, {A=1,a=1,d=1,h0=1,s=1});
    ibc:= {h(x,0)=x^2+1,h(0,t)=1,D[1](h)(0,t)=0,h(1,t)=2,D[1](h)(1,t)=2};
    S:= pdsolve(pde1,ibc,numeric);
    S:-plot3d(t=0..3, axes=box);

To remove all trailing 0's from a floating-point number, try this:

> RemTZ:= proc(x:: float)
       parse(StringTools[RegSplit]("0+$",convert(x,string)))
    end proc;

For example:

> RemTZ(123.456000);

            123.456

Just to expand a bit on Jonny's suggestion, you might look at the GetEquation Maplet application: ?Maplets,Examples,GetEquation.  See also ?examples,GetEquationMaplet

It looks like define doesn't work for partial derivatives.  You could try `diff/f`.

> `diff/f`:= proc(a,b,c) diff(a,c)*f1x(a,b) + diff(b,c)*f1y(a,b) end proc;

cos(x/n) = 0 for x = n*k*Pi/2 where k is an odd integer.  If x = m*Pi/2 for integer m (1<= m <= 2011), your function changes sign at x if and only if there is an odd number of divisors n of m for which m/n is odd, i.e. an odd number of odd divisors of m.  So:

> f:= n -> type(nops(select(type,divisors(n),odd)),odd);
   nops(select(f, [$1..2011]));

75

Somehow a complex number has appeared where only reals ought to be.  Can you upload the .bmp file, or at least give us a link to it?