with(VectorCalculus):

SurfaceInt(f,  [t1,t2,t3]=Surface(<s,t,1-s-t>,s=0..1,t=0..1-s));

I don't think Maple can find a closed form for this integral, even with appropriate assumptions on a1,...,b3. 

You might try this:

> showvals := proc(expr::uneval)
     local V, v;
     V:= indets(expr,name);
     subs(map(t -> (t = convert(eval(t),name)), V), expr)
    end proc:

> A := 10: B:= 25: C := 4:

(You really don't mean D = 6, because D is protected)

> showvals((A+B*sqrt(C))/(A+C));

(`10`+`25`*sqrt(`4`))/(`10`+`4`)

To see a phase plane with arrows, you want to convert your DE to a system of first-order DE's.  Then use phaseportrait or DEplot from the DEtools package.

> with(DEtools):
   sys := {diff(x(t),t) = xdot(t), diff(xdot(t),t) = -x(t)}:
   DEplot(sys, {x(t), xdot(t)}, t = 0 .. 7,   {[x(0)=2, xdot(0)=0]},
    x = -3 .. 3, xdot = -3 .. 3, stepsize = 0.1, linecolour=blue);

For a function of one variable, you might try allsolve in my Maple Advisor Database, www.math.ubc.ca/~israel/advisor.

For example, to find the roots of x = tan(x) for -10 < x < 10:

>  allsolve(x*cos(x) = sin(x), x=-10..10);

-7.725251837, -4.493409458, .8587267336e-4, 4.493409458, 7.725251837

Actually, Maple is smart enough to do the substitutions automatically. 

> unapply(y[10], seq(x[i],i=1..10));

(x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9, x_10) -> x_10+x_1^2+x_2*x_5+x_3*x_8

So, what are your answers?

Also try

convert(expr, FormalPowerSeries, x);

This works for some cases that convert(..., Sum) doesn't handle.  For example:

> convert(sqrt((1-sqrt(1-x))/x),Sum, x);

((1-(1-x)^(1/2))/x)^(1/2)

> convert(%, FormalPowerSeries, x);

Sum(1/2*2^(1/2)*(4*k)!*16^(-k)/(2*k)!^2/(2*k+1)*x^k,k = 0 .. infinity)

> Q := Diff(u*v, t);
   with(Student[Calculus1]):
   Rule[product](Q);

Diff(u*v,x) = Diff(u,x)*v+u*Diff(v,x)

One thing I can say:

Int(1/x/(exp(x)-1),x=r..infinity) = sum(int(1/x*exp(-n*x),x=r..infinity),n=1..infinity);

Int(1/(x*(exp(x)-1)),x = r .. infinity) = sum(Ei(1,n*r),n = 1 .. infinity)

The other terms present no problem, and it all boils down to showing that

sum(Ei(1,n*r),n=1..infinity)=r^(-1)+(1/2*gamma+1/2*ln(r)-1/2*ln(2*Pi))+o(1)

 as r -> 0+.  I haven't been able to do that...

See Quantile in the Statistics package.

> with(Statistics);  Digits:= 20;
   Quantile(FRatio(4,3),0.95);

             9.1171822532464231791

>  Quantile(StudentT(5),0.9);

           1.4758840488244810784

There are no real solutions.  Consider:

eqs:= { c1*r - M1 = c3, c1*r - M2 = c4, c2*r - M3 = c5,  r^2 = M1^2 + M2^2 + M3^2};

S := solve(eqs, {M1, M2, M3, r});

S := {M1 = c1*RootOf((-1+2*c1^2+c2^2)*_Z^2+(-2*c1*c3-2*c1*c4-2*c2*c5)*_Z+c5^2+c4^2+c3^2)-c3, M2 = c1*RootOf((-1+2*c1^2+c2^2)*_Z^2+(-2*c1*c3-2*c1*c4-2*c2*c5)*_Z+c5^2+c4^2+c3^2)-c4, M3 = c2*RootOf((-1+2*c1^2+c2^2)*_Z^2+(-2*c1*c3-2*c1*c4-2*c2*c5)*_Z+c5^2+c4^2+c3^2)-c5, r = RootOf((-1+2*c1^2+c2^2)*_Z^2+(-2*c1*c3-2*c1*c4-2*c2*c5)*_Z+c5^2+c4^2+c3^2)}

allvalues(eval(%,{c1 = 2.042662651, c2 = .8093568996, c3 = 60.36864083, c4 = 158.4676822, c5 = 101.8720814}));

{M1 = 74.81979207+46.20527768*I, M2 = -23.2792493+46.20527768*I, M3 = -48.30685327+18.30775154*I, r = 66.18245691+22.62012166*I}, {M1 = 74.81979207-46.20527768*I, M2 = -23.2792493-46.20527768*I, M3 = -48.30685327-18.30775154*I, r = 66.18245691-22.62012166*I}

In order for there to be real solutions, the following quantity has to be nonnegative:

> discrim(op(eval(r,S)),_Z);

8*c1^2*c3*c4-8*c5^2*c1^2-4*c1^2*c3^2-4*c1^2 *c4^2+8*c2*c5*c1*c3+8*c2*c5*c1*c4-4*c3^2*c2^2-4*c4^2 *c2^2+4*c3^2+4*c4^2+4*c5^2

The command is dsolve.  The derivatives can be entered as either diff(x(t),t$2) and diff(x(t),t) or (D@@2)(x)(t) and D(x)(t), or (in 2D math input) x''(t) and x'(t) (that's two single quotes, not a double quote).

I might do something like this: f(x) then VerticalLine from the Operators palette, then fraczerothick from the Layout palette, and replace the A and B by the bounds.  The result looks like this:

It doesn't work as Maple input, but it should be OK as embedded math in text.

You were plotting sqrt(1-x^2 - y), but I suspect you meant sqrt(1-x^2 - y^2).

> plot3d(sqrt(1-x^2-y^2),x=-1..1,y=-sqrt(1-x^2) .. sqrt(1-x^2));

Actually there's a secondary issue: because of roundoff error, (1-x^2-y^2) will not be exactly 0 when y = (+/-) sqrt(1-x^2), and it may be slightly negative, which will spoil the plot for those points.  So it's better to do this:

> plot3d(sqrt(max(0,1-x^2-y^2)),x=-1..1,
    y=-sqrt(1-x^2) .. sqrt(1-x^2), scaling=constrained);

This is not a worksheet (.mw) file, it's an html file containing images from the worksheet.  Please press the green up-arrow button in the editor and upload the actual worksheet.