You could use cylindrical coordinates.

> plot3d([r,theta,r^2], r = 0 .. 1, theta = 0 .. 2*Pi, 
     coords = cylindrical);

Or if you want to stick to cartesian coordinates

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

Actually there is a "pre-packaged" command of sorts.

To have Maple determine when a (scalar) expression is real, you might try

> solve(expr > -infinity);

You could try a Laplace transform in the y variable.

> heat := (1-z^2)*(diff(T(y, z), y)) = 
         (diff(z*(diff(T(y, z), z)), z))/z;

> inttrans[laplace](heat, y, s);

T(0,z)*(-1+z^2)+s*laplace(T(y,z),y,s)-z^2*s*laplace(T(y,z),y,s) = 1/z*diff(laplace(T(y,z),y,s),z)+diff(laplace(T(y,z),y,s),`$`(z,2))

One of your boundary conditions is T(0,z) = 0.

ode:= eval(%,{T(0,z)=0,laplace(T(y,z),y,s)=L(z)});

The other two will correspond to

bcs:= L(1) = inttrans[laplace](1,y,s), D(L)(0) = 0;

dsolve({ode, bcs});

L(z) = 1/exp(-1/2*I*s^(1/2))/KummerM(1/2-1/4*I*s^(1/2),1,s^(1/2)*I)/s*exp(-1/2*I*s^(1/2)*z^2)*KummerM(1/2-1/4*I*s^(1/2),1,s^(1/2)*z^2*I)

So your T(y,z) should be the inverse Laplace transform of this.  Unfortunately inttrans[invlaplace] can't do anything with it.

By saying you want it "exactly", are you implying that x is in fact the square of a rational number?  If so, sqrt(x) should find it for you.  For example:

>sqrt(72828384888321797968550560785066197301309605172436595928763160641878955647376);

  269867346836036830609348673569683065876

Here's one way to make such a grid:

> with(plots):
Grid:= seq(plot(j,theta=0..2*Pi,coords=polar,colour=gray),j=1..5),
 seq(plot([r,Pi*j/6,r=0..5], coords=polar, colour=gray), j=0..11),
 seq(textplot([j,-0.3,j]), j=1..5),
 seq(textplot([5.5, Pi*j/6, j*30], coords=polar), j=0..11):

Then

display([ your_plot, Grid], axes  = none);

Note btw that the third boundary condition is wrong: it evaluates as 0 = 0 because T(y,0) contains no z.  You mean D[2](T)(y,0) = 0.

The solution involving the Whittaker functions do not in general satisfy any of the boundary conditions.  I don't know why Maple 9.52 returned that solution.  Note, however, that without the numeric option, pdsolve does not solve boundary-value problems.  With the numeric option, it seems to work:

> heat := (1-z^2)*(diff(T(y, z), y)) = 
         (diff(z*(diff(T(y, z), z)), z))/z;
  cond := T(0, z) = 0, T(y, 1) = 1, D[2](T)(y, 0) = 0;
  Sol := pdsolve(heat,{cond},numeric);
  Sol:-plot3d(y=0..1,z=0..1,axes=box);

This is homework that I assigned.  Please don't do it for this student.

> Groebner[Basis](sys, plex(a0,a1,a2,c1,c2,m,n,k));

[k, n, c2, c1-m, a2, a1, a0-1]

So: no solutions unless a0 = 1, a1 = a2 = 0, c2 = 0, and then the only solution is m=c1, n=0, k=0.

> plot(R, 0..2*Pi, coords=polar, numpoints=N+1, adaptive=false,  
            colour=blue, filled=true);

and use display to combine this with your circle.

arcsin(x) is real for -1 <= x <= 1, but Maple is dealing with the complex version of arcsin, which actually defined everywhere in the complex plane.  In particular, for x > 1,  arcsin(x) = pi/2 - i ln(x + sqrt(x^2-1)), and therefore the limit of arcsin(x)/x  as x -> +infinity is indeed 0.

with(plots):
tplots:= map(t -> textplot([op(t+[0,0.5]),sprintf("%a",t)]), 
  data);
display(op(tplots),
  plot(data, style = point, color = black, symbol = BOX), 
  plot(cubfit, x = 0 .. 27, thickness = 2)); 
    

inequal is for systems of linear inequalities.  sin(y) is not linear.

implicitplot(sin(y),x=0..10,y=-10..10,filledregions=true,coloring=[white,green]);

To switch axes:

> plottools[transform]((x,y) -> [y,x])( your_plot );

To plot c[hat] at various times, if c[hat] depends on t and z, and z goes from a to b:

> times := [0, 1, 2, 3 ];  or whatever
    plot([seq(c[hat], t = times)], z = a .. b);

One way is like this.  It constructs a pseudo-random Vector where each entry has probability 0.17 of being 1, 0.83 of being 0.

> X := Statistics[Sample](Bernoulli(0.17), 1000);

You're thinking of densityplot, not DensityPlot.  Maple is case-sensitive.

The command you want is listdensityplot in the plots package.