In general for a nonlinear differential equation finding the radius of convergence is difficult, and a closed form for the coefficients is impossible.  For a linear differential equation, the radius of convergence is (at least, and usually exactly) the distance from the origin to the nearest singularity in the complex plane of the coefficients of the equation. 

In general for a nonlinear differential equation finding the radius of convergence is difficult, and a closed form for the coefficients is impossible.  For a linear differential equation, the radius of convergence is (at least, and usually exactly) the distance from the origin to the nearest singularity in the complex plane of the coefficients of the equation. 

Sorry, of course I meant 3/2*b = b + 1/(1/b + 1/b)

In this case the red and blue regions wouldn't overlap, but what would you want to get if they did?

Here's one solution (when they don't overlap).

with(plots):
F1 := min(1-(x-y), x+y);
F2 := min(-x, -y);
P1:= select(has, implicitplot(F1,x=-3..3,y=-3..3,gridrefine=3,
          filledregions=true,coloring=[white,red]),
      COLOUR(RGB, 1.00000000, 0., 0.)):
P2:= select(has, implicitplot(F2,x=-3..3,y=-3..3,gridrefine=3,
          filledregions=true,coloring=[white,blue]),
      COLOUR(RGB, 0., 0., 1.00000000)):
 display(P1,P2);

In this case the red and blue regions wouldn't overlap, but what would you want to get if they did?

Here's one solution (when they don't overlap).

with(plots):
F1 := min(1-(x-y), x+y);
F2 := min(-x, -y);
P1:= select(has, implicitplot(F1,x=-3..3,y=-3..3,gridrefine=3,
          filledregions=true,coloring=[white,red]),
      COLOUR(RGB, 1.00000000, 0., 0.)):
P2:= select(has, implicitplot(F2,x=-3..3,y=-3..3,gridrefine=3,
          filledregions=true,coloring=[white,blue]),
      COLOUR(RGB, 0., 0., 1.00000000)):
 display(P1,P2);

Interesting... 

 convert(Ei(1,x), GAMMA);

GAMMA(0,x)

(which is correct).

 convert(Ei(1,x), FormalPowerSeries, x);

-gamma-ln(x)+Sum((-1)^k/(k+1)!*x^(k+1)/(k+1),k = 0 .. infinity)

(also correct, I think).

Why is the <maple> tag not working for me?

Here's a way to get Maple to discover a solution.  Of course it won't be a very nice solution...

First, here's a procedure to generate random expressions of maximum depth n using a,b,1,+,-,and reciprocal.
 

  randt:= proc(n) 
      local r,q;
      if n = 0 then 
        r := rand(1..3)();
        if r = 1 then a elif r = 2 then b else 1 fi
      else
       do 
        r := rand(1..7)();
        if r = 1 then return a
        elif r = 2 then return b
        elif r = 3 then return 1
        elif r = 4 then q:= randt(n-1) + randt(n-1);
        elif r = 5 then q:= randt(n-1) - randt(n-1);
        else return 1/randt(n-1)
        fi;
        if normal(q) <> 0 then return q fi
       end do;
      fi
 end proc;

Gather results that are equivalent to polynomials of degree >= 2.  Make sure you get one that has a term in a*b.

count := 0; going:= true;
while going do
     T := randt(6);
     N := normal(T);
     if type(N, polynom) and degree(N) >= 2 then
       count:= count+1;
       Exprs[count]:= T;
       NExprs[count]:= expand(N);
       if coeff(coeff(%,a),b) <> 0 then going:= false fi;
       print(T,N);
     end if
 end do:

Now see if a linear combination of these results will do it.

Q:=collect(add(x[i]*NExprs[i],i=1..24)+x[0]+x[a]*a+x[b]*b,
     [a,b],distributed);
S:= solve({coeffs(Q - a*b, [a,b])});
arbitrary := map(t -> (t=0), indets(map(rhs,S)));
S1 := eval(S, arbitrary);
eval(add(x[i]*Exprs[i],i=1..24)+x[0]+x[a]*a+x[b]*b, 
     S1 union arbitrary);

3/2*b+3/2*a+2/(1/(b-4)-1/b)+1/(2*(1/(b-a)-1/(-a+b-1)))+1/(1/(a-2)-1/a)

Note that e.g. 3/2*b = b + 1/(b+b)

If you're using Windows, one of the Maple 11 items you can find in the Start Menu or on your desktop should be called "Classic Worksheet Maple 11".  If you're using Unix/Linux, the command is

maple -cw

On the Macintosh there is no Classic.

If you're using Windows, one of the Maple 11 items you can find in the Start Menu or on your desktop should be called "Classic Worksheet Maple 11".  If you're using Unix/Linux, the command is

maple -cw

On the Macintosh there is no Classic.

You could write the solution as

1/(2*(1/(a+b)-1/(a+b+1)))-1/(2*(1/a-1/(a+1)))-1/(2*(1/b-1/(b+1)))

where 2*x means x+x.

I didn't need to figure it out: I first entered it from the palette and then used Convert to, 1D Math Input.

I didn't need to figure it out: I first entered it from the palette and then used Convert to, 1D Math Input.

Here's a nice example, I think.

> de := diff(y(x),x$2) = 
      -2.645751311/(7.*x^4-5.291502622*x^2+1.)*y(x);

> dsolve(de);

y(x) = _C1*HeunG((-129181180998500967-48825897*I*341781279^ (1/2))/(-129181180998500967+48825897*I*341781279^(1/2)), -2645751311/195303588*I/(2645751311+341781279^(1/2)*I) *(-48825897+377964473*I*341781279^(1/2))*341781279^(1/2) /(-2645751311+341781279^(1/2)*I),-1/2,0,1/2,0, -7000000000/(-2645751311+341781279^(1/2)*I)*x^2) +_C2*x*HeunG((-129181180998500967-48825897*I *341781279^(1/2))/(-129181180998500967+48825897*I *341781279^(1/2)),-2645751311/195303588*I/ (2645751311+341781279^(1/2)*I)*(-48825897+377964473*I* 341781279^(1/2))*341781279^(1/2)/(-2645751311+ 341781279^(1/2)*I),0,1/2,3/2,0,-7000000000/ (-2645751311+341781279^(1/2)*I)*x^2)

> dsolve(convert(de, rational));

y(x) = _C1*HeunG(16129/16128,-32257/129024,-1/2,0,1/2,0, 127/48*x^2)+_C2*HeunG(16129/16128,-32257/129024,0, 1/2,3/2,0,127/48*x^2)*x

> dsolve(identify(de));

y(x) = _C1*(7^(1/2)*x^2-1)^(1/2)+_C2*(7^(1/2)*x^2-1)^(1/2) *ln((7^(1/4)*x+1)/(7^(1/4)*x-1))

The general idea is that you get a sequence of solutions satisfying these (homogeneous) boundary conditions, and the general solution will be written as a
series in these basic solutions.  You won't get exponentials to be 0, but a trig function (which is the real or imaginary part of a complex exponential) can be.  In this case
exp(k1*x+k2*t) is replaced by sin(Pi*n*x)*exp(k2*t).  Thus:

> Ansatz:= {u1(x,t)=A*sin(n*Pi*x)*exp(k2*t), u2(x,t)=B*sin(n*Pi*x)*exp(k2*t),
                     u3(x,t)=C*sin(n*Pi*x)*exp(k2*t)};
  

and

> ABCeqs:= eval(eval(des, Ansatz), {exp=1,sin(n*Pi*x)=1});

Usually you'll also have initial conditions, typically U(x,0) = F(x) for 0 <= x <= 1,
and you use sine Fourier series for the components of S^(-1).F(x) (in the notation
I used in "More general approach") to get the coefficients for the basic solutions.
Just about every DE textbook that covers the heat equation will show how to do
this in the case of a single heat equation; it's essentially the same here.

Are you asking how I got the Ansatz?  As the name implies, it's an educated guess rather than a Maple command (see en.wikipedia.org/wiki/Ansatz). 
 

BTW, I'm rather surprised that pdsolve won't give any result for this system.  It would produce a (rather complicated) solution if at least two of the constants a,b,c were given values.  Try, for example,

> pdsolve(eval(des, {a=1,b=2,c=3}));