Certain automatic simplifications are applied to all output, and are very hard to avoid.  This includes taking out that factor of 1/3.  Well, in a pinch you could do this:

> (5*x+4)/Typesetting:-Typeset(3*x) = 7;

OK, there is indeed a new and serious bug in matrix multiplication in Maple 15, affecting non-square diagonal Matrices.  Consider:

> A:= Matrix(3,2,symbol=a, shape=diagonal);
   B:= Matrix(2,2,symbol=b);
   A . B;
                 [a[1, 1] b[1, 1]    a[2, 2] b[2, 2]]
                 [                                  ]
                 [a[2, 2] b[2, 1]           0       ]
                 [                                  ]
                 [a[1, 1] b[1, 2]           0       ]

It should, of course, be 

                 [a[1, 1] b[1, 1]    a[1, 1] b[1, 2]]
                 [                                  ]
                 [a[2, 2] b[2, 1]    a[2, 2] b[2, 2]]
                 [                                  ]
                 [       0                  0       ]

Even worse:

> A:= Matrix(2,3, symbol=a, shape=diagonal);
   C:= Matrix(3,3, symbol=c);
   A . C;

mserver.exe has stopped working

I have submitted an SCR.

@J F Ogilvie : could you provide more details about this bug?  In particular, what are SI and V?
For ordinary Matrices (of compatible dimensions) LinearAlgebra:-Transpose(SI) . V works fine as far as I can tell.   

See the help page ?plot,options.
The plot command computes a certain number of points on your curve.  By default it uses adaptive plotting, so that more points are computed where the function appears to be changing rapidly.  Now if you use style=point, Maple will plot the points it calculates.  On the other hand, linestyle=dot plots equally-spaced dots on the line segments joining those points.

If you want equal spacing in the x direction with the style=point option, you can turn off adaptive plotting with the option adaptive=false.  Equal spacing with respect to arc length in the point style would not be easy to arrange. 

A cute idea, but the formula you display is not the actual formula you used to produce it.

> tupper:= simplify(tupper,size) assuming x::posint, y::posint;

The next question is how to produce the appropriate n for this expression.  It would be nice to do that all within Maple, but that doesn't seem possible.  In Maple you can export a worksheet to html with the formulas exported as gif files.  But Maple's ImageTools package can't handle gif, only jpeg, tiff and bmp.  I'll have to go outside Maple to convert a gif file to one of those.  I used Paint to do this (in this case saving as jpeg).  While we're at it, I changed the 17 to 64, because the jpeg file turns out to have 64 rows. 

>   R:= ImageTools:-Read("d:/test/images/tupper_1.jpg");

The resut is a 1..64 x 1..569 x 1..3 Array of float[8].  We want to make this into an integer with 1 bits encoding positions where this is dark.

> N:= add(add(2^(64*j+i)*piecewise(R[i,j,1]>0.6 , 0, 1),i=1..64),j=1..569):

3063406762436416313461751591304266739808654845271255499869500714894499452810708430119355129597633782[...10581 digits...]2606770365668881384818361357051076329155328042582836994679242715489051568901043552066927768960499712

Now to produce the array of values using the expression tupper (or more precisely, a function corresponding to it).

> tf:= unapply(tupper,x,y);
   A:= Array(1..64,1..569,(y,x)-> tf(x,y+N),datatype=float[8],order=C_order):

This may be a bit big for listdensityplot to handle comfortably; better to use ImageTools:-Create.

> Img:= ImageTools:-Create(A);
   ImageTools:-View(Img);

Your procedure would generate uniformly-distributed random points on a sphere, but when you map this to an ellipsoid the distribution may no longer be uniform with respect to area.  Suppose you generate random points for parameters (u,v) in a rectangle (in your case [-1,1] x [0,2*Pi]) and map (u,v) to F(u,v) in 3-dimensional space.  The mapping multiplies areas by the length of the cross product of the vectors and ; call this C(u,v).  In order for the mapped points to be uniformly distributed with respect to area, you want the density of your points in (u,v) space to be proportional to C(u,v).  In your case

F:= (u,v) -> <rX*sqrt(1-u^2)*cos(v), rY*sqrt(1-u^2)*sin(v), rZ*u>;
CP:= VectorCalculus:-CrossProduct(diff~(F(u,v),u), diff~(F(u,v),v));
C:= unapply(simplify(sqrt(CP[1]^2 + CP[2]^2 + CP[3]^2)),u,v);

We can obtain our (u,v) points by a rejection sampling method: start out with uniformly distributed points, and accept or reject each one, accepting with probability proportional to C(u,v).  If 0 < rX <= rY <= rZ, the maximum of C is at u=0, v=0, where C(0,0) = rY*rZ.  So we accept (u,v) with probability C/(rY*rZ).  For example:

> rX, rY, rZ:= 0.5, 1, 1.5;
  C0:= rY*rZ;
  with(Statistics):
  RawPts:= zip(`[]`,Sample(Uniform(-1,1),10000),Sample(Uniform(0,2*Pi),10000)):
  AVals:= Sample(Uniform(0,1),10000):
  Accepted:= select(i -> (AVals[i] < C(op(RawPts[i]))/C0), [$1..1000]):
  uvPts:= RawPts[Accepted];
  Pts:= map(p -> [rX*sqrt(1-p[1]^2)*cos(p[2]),rY*sqrt(1-p[1]^2)*sin(p[2]),rZ*p[1]], uvPts);
  plots[pointplot3d](Pts,axes=box,scaling=constrained,colour=red);

I suspect that line was introduced to correct a different bug.  The mod function behaves a bit strangely with modulus of 1:  t mod 1 returns t, although c*t mod 1 (for any integer c other than 1) returns 0.  As a result, chrem([t, 1+t^2], [1,3]) returns 3*t + 1 + t^2 before Maple 13 and 1 + t^2 after Maple 13.

I suspect that line was introduced to correct a different bug.  The mod function behaves a bit strangely with modulus of 1:  t mod 1 returns t, although c*t mod 1 (for any integer c other than 1) returns 0.  As a result, chrem([t, 1+t^2], [1,3]) returns 3*t + 1 + t^2 before Maple 13 and 1 + t^2 after Maple 13.

To get the inert version, you can type in

Diff(x^2*sin(x),x)

To get the inert version, you can type in

Diff(x^2*sin(x),x)

Indeed.  In general, if (as Maple does) you use the principal branches of powers,
(a^(x+2))^(1/3) = exp((x+2)*ln(a) + 2*Pi*I*m/3) where m = floor(1/2 - Im((x+2)*ln(a))/(2*Pi)),
and similarly
(a^(x-5))^(1/2) = exp((x-5)*ln(a) + 2*Pi*I*n/2) where n = floor(1/2 - Im((x-5)*ln(a))/(2*Pi))).
The equation is true if (19-x)/6 * ln(a) = 2*Pi*I*(n/2 - m/3 + k) for some integer k.
If a = -1, it seems there are no solutions. 

Indeed.  In general, if (as Maple does) you use the principal branches of powers,
(a^(x+2))^(1/3) = exp((x+2)*ln(a) + 2*Pi*I*m/3) where m = floor(1/2 - Im((x+2)*ln(a))/(2*Pi)),
and similarly
(a^(x-5))^(1/2) = exp((x-5)*ln(a) + 2*Pi*I*n/2) where n = floor(1/2 - Im((x-5)*ln(a))/(2*Pi))).
The equation is true if (19-x)/6 * ln(a) = 2*Pi*I*(n/2 - m/3 + k) for some integer k.
If a = -1, it seems there are no solutions. 

In both of your equations, one side depends only on t and the other side depends only on r.  The only way that can be true is if both sides are constant.  So both u(r) and v(t) are constants, and it's not hard to see that these constants must both be 0.  Of course this doesn't satisfy your boundary conditions.

@Christopher2222 :

So does D[1](T)(2,t)=0 mean there that the rate of heat change there is zero?

No, it's the derivative with respect to x rather than t, so it means that the temperature gradient is 0.  And since the temperature gradient is proportional to the heat flux, that means that heat is not flowing across the boundary.  

One thing I find is that if I use a 2d plot at a specific time say t=100000, Maple plods through calculations up to that point before it presents me with the graph.  Shouldn't that be faster, I mean nearly instantaneous?

Maple is using numerical methods here, rather than a symbolic solution.  In order to find the solution at t=100000, it must start with the initial conditions at t=0 and see how things change as time progresses, one time step at a time.

And how can I get the value at a specific time and point on the rod, say at x=2 and t=10800 (3 hours)?  Skin burns at around 130 F or about 50 C (323 K) and after 3 hours, from the graph, it is roughly just below 330K probably still too hot to touch the end of the rod.  I couldn't make sense of the help pages to determine a value at a specific point in time.

That can be done using value.  Thus:

> v:= sol:-value():
   v(0.2,3); 
  

@Christopher2222 :

So does D[1](T)(2,t)=0 mean there that the rate of heat change there is zero?

No, it's the derivative with respect to x rather than t, so it means that the temperature gradient is 0.  And since the temperature gradient is proportional to the heat flux, that means that heat is not flowing across the boundary.  

One thing I find is that if I use a 2d plot at a specific time say t=100000, Maple plods through calculations up to that point before it presents me with the graph.  Shouldn't that be faster, I mean nearly instantaneous?

Maple is using numerical methods here, rather than a symbolic solution.  In order to find the solution at t=100000, it must start with the initial conditions at t=0 and see how things change as time progresses, one time step at a time.

And how can I get the value at a specific time and point on the rod, say at x=2 and t=10800 (3 hours)?  Skin burns at around 130 F or about 50 C (323 K) and after 3 hours, from the graph, it is roughly just below 330K probably still too hot to touch the end of the rod.  I couldn't make sense of the help pages to determine a value at a specific point in time.

That can be done using value.  Thus:

> v:= sol:-value():
   v(0.2,3);