Oops, that was the identity that comes out of of the integral of sqrt(1+x^(-2/3)).
This one will be a bit different, but I think the same principles should apply.

> J := int(1/3*(9+2^(2/3)/x^(2/3))^(1/2),x = 0 .. U) assuming U > 0;

>  factor(27*J + 2);

(9*U^(2/3)+2^(2/3))^(3/2)

> convert(eval(%, U = t^(-1/2)), FormalPowerSeries, t);

Sum(-27*(-1)^k*pochhammer(-1/6,k)*pochhammer(1/6,k)/k!/(3*k)!*(2*k)!/(2*k-1)*27^(-k)*t^(k-1/2),k = 0 .. infinity)+Sum(9/2*2^(2/3)*(-1)^k*pochhammer(-1/6,k)*(2*k)!*pochhammer(1/6,k)*27^(-k)/k!/(3*k+1)!*t^(-1/6+k),k = 0 .. infinity)+Sum(1/2*2^(1/3)*(-1)^k*pochhammer(1/6,k)*pochhammer(5/6,k)*(2*k)!*27^(-k)/k!/(3*k+2)!*t^(1/6+k),k = 0 .. infinity)

> map(convert, %, hypergeom);

27/t^(1/2)*hypergeom([-1/2, -1/6, 1/6],[1/3, 2/3],-4/729*t)+9/2*2^(2/3)/t^(1/6)*hypergeom([-1/6, 1/6, 1/2],[2/3, 4/3],-4/729*t)+1/4*2^(1/3)*t^(1/6)*hypergeom([1/6, 1/2, 5/6],[4/3, 5/3],-4/729*t)

> simplify(eval(J, U = 2) = eval((%-2)/27, t = 1/4));

-2/27+20/27*10^(1/2) = 253/108*hypergeom([-1/2, -1/6, 1/6],[4/3, 5/3],-1/729)-911/6298560*hypergeom([1/2, 5/6, 7/6],[7/3, 8/3],-1/729)+73/9795520512*hypergeom([3/2, 11/6, 13/6],[10/3, 11/3],-1/729)-2/27

There are three main parts of Maple: the kernel, the library, and the user interface.
See e.g. page 396 of the Maple 12 Advanced Programming Guide www.maplesoft.com/view.aspx for a description of what they do.

There are three user interfaces: command-line, Classic and Standard.  Standard is the one you're almost certainly using if you haven't made an effort to use Classic or command-line.  Yes, the appearance of the worksheet is part of what the user interface handles.

The identity seems to be Eq = 0 where

> Eq := 3/8*1/u^(1/3)*hypergeom([1/6, 1/2, 5/6],[4/3, 5/3],-1/(u^2))
+3/2*u^(1/3)*hypergeom([-1/6, 1/6, 1/2],[2/3, 4/3],-1/(u^2))
+u*hypergeom([-1/2, -1/6, 1/6],[1/3, 2/3],-1/(u^2))
- (u^(2/3)+1)^(3/2);

Change variables so -1/u^2 becomes t:

> eval(Eq, u = t^(-1/2));

Expand as power series in t:

> map(convert, %, FormalPowerSeries, t);

> expand(%);

                                            0

The identity seems to be Eq = 0 where

> Eq := 3/8*1/u^(1/3)*hypergeom([1/6, 1/2, 5/6],[4/3, 5/3],-1/(u^2))
+3/2*u^(1/3)*hypergeom([-1/6, 1/6, 1/2],[2/3, 4/3],-1/(u^2))
+u*hypergeom([-1/2, -1/6, 1/6],[1/3, 2/3],-1/(u^2))
- (u^(2/3)+1)^(3/2);

Change variables so -1/u^2 becomes t:

> eval(Eq, u = t^(-1/2));

Expand as power series in t:

> map(convert, %, FormalPowerSeries, t);

> expand(%);

                                            0

The fact that integrals such as this produce a complicated result using hypergeometric functions when a simple elementary antiderivative is known is, I think, a bug.  For example:

> int((1+1/x^(2/3))^(1/2), x=0..2);

-1/48*2^(2/3)/Pi^(5/2)*3^(1/2)*(-3*Pi^(5/2)*3^(1/2)*hypergeom([1/6, 1/2, 5/6],[4/3, 5/3],-1/4)+8*Pi^(5/2)*3^(1/2)*2^(1/3)-12*Pi^(5/2)*3^(1/2)*2^(2/3)*hypergeom([-1/6, 1/6, 1/2],[2/3, 4/3],-1/4)-16*Pi^(5/2)*3^(1/2)*2^(1/3)*hypergeom([-1/2, -1/6, 1/6],[1/3, 2/3],-1/4))

> eval(int((1+1/x^(2/3))^(1/2), x=0..r), r=2) assuming r > 0;

-1+(2^(2/3)+1)^(1/2)*2^(2/3)+(2^(2/3)+1)^(1/2)

The fact that integrals such as this produce a complicated result using hypergeometric functions when a simple elementary antiderivative is known is, I think, a bug.  For example:

> int((1+1/x^(2/3))^(1/2), x=0..2);

-1/48*2^(2/3)/Pi^(5/2)*3^(1/2)*(-3*Pi^(5/2)*3^(1/2)*hypergeom([1/6, 1/2, 5/6],[4/3, 5/3],-1/4)+8*Pi^(5/2)*3^(1/2)*2^(1/3)-12*Pi^(5/2)*3^(1/2)*2^(2/3)*hypergeom([-1/6, 1/6, 1/2],[2/3, 4/3],-1/4)-16*Pi^(5/2)*3^(1/2)*2^(1/3)*hypergeom([-1/2, -1/6, 1/6],[1/3, 2/3],-1/4))

> eval(int((1+1/x^(2/3))^(1/2), x=0..r), r=2) assuming r > 0;

-1+(2^(2/3)+1)^(1/2)*2^(2/3)+(2^(2/3)+1)^(1/2)

To some extent it may depend on what you're trying to do with Maple.  But for ordinary
purposes, I think 512 MB should be plenty.

The simplest way to define your function is

> f := x -> piecewise(x>=r, sqrt(x-r), -sqrt(r-x));

The Newton procedure for one iteration:

> Newt:= x -> evalf(x - f(x)/D(f)(x));

Now the way I like to visualize this is by drawing a vertical line from [x[i], 0] to [x[i], f(x[i])] and then the tangent line from there to [x[i+1],0].  The following procedure will do this.

> NewtStep:= x -> plot([[x,0],[x,f(x)],[Newt(x), 0]], colour=blue);

In your particular case, it turns out that the Newton function is an involution: Newt(Newt(x)) = x.  Here's a plot:

> r:= 1; 
  plots[display]([plot(f(x),x=-1..3), NewtStep(2.5),
    NewtStep(Newt(2.5))]);

The simplest way to define your function is

> f := x -> piecewise(x>=r, sqrt(x-r), -sqrt(r-x));

The Newton procedure for one iteration:

> Newt:= x -> evalf(x - f(x)/D(f)(x));

Now the way I like to visualize this is by drawing a vertical line from [x[i], 0] to [x[i], f(x[i])] and then the tangent line from there to [x[i+1],0].  The following procedure will do this.

> NewtStep:= x -> plot([[x,0],[x,f(x)],[Newt(x), 0]], colour=blue);

In your particular case, it turns out that the Newton function is an involution: Newt(Newt(x)) = x.  Here's a plot:

> r:= 1; 
  plots[display]([plot(f(x),x=-1..3), NewtStep(2.5),
    NewtStep(Newt(2.5))]);

I frequently run Maple 12 on an old computer with 256 MB of RAM (PIII/601 Mhz,
Windows XP).  Yes, it is slow, and if I try to have too many windows open at once it
will complain about running out of virtual memory, but it does work.

Specifically, the problem is this:

> int(BesselK(0,r)*r, r = 0 .. b);

-BesselK(1,b)*b+1/2*b*BesselK(0,b)*BesselI(1,b)+1/2*b*BesselK(1,b)*BesselI(0,b)

A correct answer would be 1/2 + this, although a simpler correct answer would be 1 - b*BesselK(1,b).  It's rather strange because the following is correct:

> int(BesselK(0,r)*r, r = a .. b);

BesselK(1,a)*a-BesselK(1,b)*b

> limit(%, a=0);

1-BesselK(1,b)*b

I suggest you check your theory again.  Your C2 was

GAMMA(c, t/b)/GAMMA(c)

Recall that GAMMA(a, z) = int(exp(-s)*s^(a-1), s = z .. infinity)

At t=0, C2 would be 1, and as t -> infinity it would go to 0.  So this is not a CDF.  1-C2 is a CDF.

I suggest you check your theory again.  Your C2 was

GAMMA(c, t/b)/GAMMA(c)

Recall that GAMMA(a, z) = int(exp(-s)*s^(a-1), s = z .. infinity)

At t=0, C2 would be 1, and as t -> infinity it would go to 0.  So this is not a CDF.  1-C2 is a CDF.

One might quibble that u*Heaviside(u) and v*Heaviside(v) are not really solutions, because they are not differentiable everywhere wrt u and v respectively when those variables are 0.  But I'll ignore that.  You want a solution of this form:

> Form := F(u*Heaviside(u), v*Heaviside(v));

Now you should be able to use pdetest as follows:

> pdetest(f(u,v) = Form, eqn);

or pdsolve with a hint:

> pdsolve(eqn, HINT = Form);

but these both seem to run into the same bug in Maple:

Error, (in limit/sumprod) invalid input: has expects 2 arguments, but received 3
 

Looking at lasterror, it seems that the bug may actually be in simplify:

> simplify(D[1,2](F)(g(u),g(v))*u*Dirac(u));

Error, (in limit/sumprod) invalid input: has expects 2 arguments, but received 3
 

I'll send in an SCR for this.

Anyway, let's plug your form into the equation and see what we get.

> eval(eqn, f(u,v)=Form);

D[1,2](F)(u*Heaviside(u),v*Heaviside(v))*(Heaviside(v)+v*Dirac(v))*(Heaviside(u)+u*Dirac(u)) = 0

Now the last two factors simplify to 0 when u < 0 or v < 0, and we're ignoring differentiability concerns when u = 0 or v = 0, so the main question is what happens when u > 0 and v > 0.

> % assuming u > 0, v > 0;

D[1,2](F)(u,v) = 0

And to find a solution of that PDE (which of course happens to be the same as our eqn):

> pdsolve(%);

F(u,v) = _F2(u)+_F1(v)

So there you have your solution:

f(u,v) = _F2(u*Heaviside(u)) + _F1(v*Heaviside(v))

One might quibble that u*Heaviside(u) and v*Heaviside(v) are not really solutions, because they are not differentiable everywhere wrt u and v respectively when those variables are 0.  But I'll ignore that.  You want a solution of this form:

> Form := F(u*Heaviside(u), v*Heaviside(v));

Now you should be able to use pdetest as follows:

> pdetest(f(u,v) = Form, eqn);

or pdsolve with a hint:

> pdsolve(eqn, HINT = Form);

but these both seem to run into the same bug in Maple:

Error, (in limit/sumprod) invalid input: has expects 2 arguments, but received 3
 

Looking at lasterror, it seems that the bug may actually be in simplify:

> simplify(D[1,2](F)(g(u),g(v))*u*Dirac(u));

Error, (in limit/sumprod) invalid input: has expects 2 arguments, but received 3
 

I'll send in an SCR for this.

Anyway, let's plug your form into the equation and see what we get.

> eval(eqn, f(u,v)=Form);

D[1,2](F)(u*Heaviside(u),v*Heaviside(v))*(Heaviside(v)+v*Dirac(v))*(Heaviside(u)+u*Dirac(u)) = 0

Now the last two factors simplify to 0 when u < 0 or v < 0, and we're ignoring differentiability concerns when u = 0 or v = 0, so the main question is what happens when u > 0 and v > 0.

> % assuming u > 0, v > 0;

D[1,2](F)(u,v) = 0

And to find a solution of that PDE (which of course happens to be the same as our eqn):

> pdsolve(%);

F(u,v) = _F2(u)+_F1(v)

So there you have your solution:

f(u,v) = _F2(u*Heaviside(u)) + _F1(v*Heaviside(v))