Errors when thickness is reduced solving 2D cylindrical heat conduction equation

I am trying to solve the heat conduction equation in cylindrical coordinates with a defined heat generation within the volume (a thin disk, with different heat flux BC top and bottom). The problem is asimuthally symmetric so I only have derivatives in "r" and "z". I am trying to use the method discussed on this site for solving inhomogeneous PDEs. I was able to get a solution when using a thickness equal to the radius (zt=1), but when I reduce the thickness (zt=0.2) I get the following errors

NDSolve::mxst: Maximum number of 10000 steps reached at the point t == 1.5022438114283172. >>

NDSolve::eerr: Warning: Scaled local spatial error estimate of 11.936241527250862 at t = 1.5022438114283172 in the direction of independent variable r is much greater than prescribed error tolerance. Grid spacing with 15 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or you may want to specify a smaller grid spacing using the MaxStepSize or MinPoints method options. >>

How can I fix this (I am using Mathematica Ver 8)?

Code:

qdot = 10; (*heat input*)
k = 1;  (*thermal conductivity*)
ti = 20; (*iniital temp*)
hb = 1; (*convection coeff back*)
hf = .01; (*convection coeff front*)
tb = 20; (*coolant temp back*)
tf = 20; (*coolant temp front*)
zt = .2; (*thickness*)

Clear[y];
f1[r_ /; r > 0] := 1;
f1[r_ /; r == 0] := 2;
f2[r_ /; r > 0] := 1/r;
f2[r_ /; r == 0] := 0;
qdot/k + Derivative[0, 0, 2][y][t, r, z] +
f2[r]*Derivative[0, 1, 0][y][t, r, z] +
f1[r]*Derivative[0, 2, 0][y][t, r, z] ==
Derivative[1, 0, 0][y][t, r, z]

eqns = {qdot/k + Derivative[0, 0, 2][y][t, r, z] +
f2[r]*Derivative[0, 1, 0][y][t, r, z] +
f1[r]*Derivative[0, 2, 0][y][t, r, z] ==
Derivative[1, 0, 0][y][t, r, z],
Derivative[0, 1, 0][y][t, r0, z] == 0,
Derivative[0, 0, 1][y][t, r, zt] == -.5,
Derivative[0, 1, 0][y][t, 0, z] == 0,
Derivative[0, 0, 1][y][t, r, 0] == +1,
y[0, r, z] == z - z^2/(2*zt)*1.5};

y[t_, r_, z_] = y[t, r, z] /.
First[NDSolve[eqns,
y[t, r, z], {t, 0, 1.6}, {r, 0, r0}, {z, 0, zt},
Method -> {"MethodOfLines", Method -> "StiffnessSwitching",
"DifferentiateBoundaryConditions" -> {True,
"ScaleFactor" -> 0}}]];

• You may be able to obtain a solution for zt = 1, but it does not look like a credible solution to me. Clearly, there are problems at r = z = 0`, and perhaps elsewhere. – bbgodfrey Feb 19 '15 at 3:00