@tomleslie There is no requirement of continuity in boundary conditions.  Discontinuous boundary conditions are quite common, and in fact absolutely necessary, in applications of PDEs.  Discontinuous boundary conditions, as well as discontinuous solutions(!) are interpreted in functional analytic setting through what are called generalized derivatives.

Indeed, the boundary value problem posed in this thread may be solved exactly.  The solution is
u(x,t)=Heaviside(t - x/4);
Maple's numeric solver gets a reasonable approximation to this if spacestep is sufficiently fine, as in
sol:= pdsolve(eq1,cond, numeric,time=t,range=0..10,spacestep=0.01);

Plot the exact and approximate solutions together to see how they compare.

And of course u(x,t) is discontinuous at (0,0) so we don't ask of its its value there.

@nm I trust that you see that the r derivative of u at r=0 should be zero by symmetry. Now look at the general solution (7.9.21) on page 319 in Haberman (I am looking at the 3rd edition). There the r dependence enters as a weighted sum of Bessel functions of all orders.  In particular, the sum includes the Bessel function of order 1 [which is J_1(...) in mathematical notation, BesselJ(1,...) in Maple notation].  But the derivative of J_1 at zero is not zero(!) so if you leave that term in, your solution u will have a nonzero derivative at r=0 which shouldn't be there. The 4th boundary condition that I have supplied gets rid of that term. 

You said you solved that equation by hand. Check to see what you did with the J_1 term.

Edit:

I take back what I have written regarding the extra boundary condition.  I have verified that that boundary condition is automatically satisfied (after dropping the singular Bessel functions as you have noted.)  I was thrown by mixing this up with the case where u depends on the polar angle, and in which case Bessel functions of all orders show up.  In your case only the Bessel function J_0 shows up, which does have zero derivative at the origin, so it's just fine.

The boundary conditions are incomplete.  Here I will describe how to correct them, but that doesn't help with the error message which appears to be due to a bug in Maple.

What is the domain of the unknown function u(r,z,t)?   It is
0 < r < 1,    0 < z < 1,   t > 0.

Note that for a fixed t, the domain is a rectangle in the (r,z) plane.  To have a well-posed problem, you need to specify boundary conditions all around that rectangle.  Your bc gives boundary conditions on three of that rectangle's four edges.  A boundary condition on the edge r=0 is missing.  The missing condition is D[1](u)(0,z,t)=0, which follows from the symmetry of u in r.

Thus, the correct boundary conditions are:
bc := u(r,0,t)=0, u(r,1,t)=0, u(1,z,t)=0, D[1](u)(0,z,t)=0;
Other than that, the problem is posed correctly.  Maple should be able to produce a series solution for arbitrary ic := u(r,z,0) = f(r,z).  But it fails even for the simpler ic := u(r,z,0) = 1.

@cencen_cj  I don't use the debugger so I cannot be of help here.  Perhaps someone else can offer advice.

@Carl Love What is wrong with this simple way?  After isolate(coeff1, D__1)*C1  I would do:
int(rhs(%), T[1]);
 

@Joe Riel Thanks for pointing out the SMT. Someone else mentioned that to me a couple of weeks ago in a totally different context. I may look into that when I get the time.

@Carl Love Thanks for your comments.  I see things a little more clearly now.

And regarding Whitehead and Russell's Principia, if my memory serves,, they arrive at the derivation of 1+1=2 somewhere in the middle of volume 2, which discourages me from pursuing that avenue :-)

@tomleslie The originally stated boundary conditions did not make sense. It was sort of obvious to me what the intention was and I corrected them, as you did. The boundary conditions are not the culprit.

The source of the problem is that the PDE is ill-posed. Its 4th order term has the wrong sign. I can't tell why the OP is interested in that PDE but I suspect that the sign of the 4th order term is incorrect. The coefficient should be −5 instead of +5. The numerical "solution" that appears in my post looks like a nice graph but probably it is total garbage.  It's no wonder that changing the solver's spacestep parameter drastically changes the solver's output.

The following stripped-down version of the PDE illustrates the problem's ill-posed nature.

Download ill-posed-pde.mw

Joe, that's very clever. I knew nothing about the Logic package until now. I would be interested in learning more but I am having a hard time wrapping my head around it. If you get the chance, could you please show how one may formulate the following (very) elementary problem with the Logic tools?

Find integers a and b selected from the set {2,4,5} so that b > a and a^2 + b^2 > 28.

There are two solutions: {a=2,b=5} and {a=4,b=5}, but I will be happy with either one.

I have no idea whether the Logic package is suitable for such a problem.

@vv Smoothing out the discontinuities certainly is an option.  In fact, that was the first solution that I offered when this came up in a recent Question. The problem with smoothing is that it requires an exceedingly fine mesh to achieve good accuracy as we see in your choice of timestep=0.0001. It was the desire to avoid such fine meshes that motivated me to seek an alternative approach which then led to the subject of this Post.

@vv Here is an exact solution corresponding to discontinuous conductivity and discontinuous initial data.

I have added this example to the previous worksheet.  Download updated worksheet: heat-finite-difference-v3.mw

@vv and @Bart, I have modified the original worksheet by adding a third example in which I compare the finite difference solution against the exact solution corresponding to a problem with discontinuous initial data.  This animation, extracted from that worksheet, demonstrates that the exact and numerical solutions match quite well.  The graphs of the exact solution (blue) and the finite difference solution (red) are difficult to distinguish.

I have not made an error analysis to determine the accuracy of the method.

By the way, from the geeral theory of parabolic equations we know that although the initial condition is discotiuous, the solution at any t>0 is in C^∞.

Download updated worksheet:  heat-finite-difference-v2.mw

@Carl Love and acer, thank your for all your detailed comments on solving tridiagonal systems.  It is now obvious to me that inverting a tridiagonal matrix for solving a linear system is a very bad idea.  That the system needs is solved repeatedly within a for-loop does still not justify inverting the matrix.  I stand corrected.

@vv That's excellent!  Converting the 2nd order PDE into a system of 1st order PDEs does the trick.  In fact, that's what I did in deriving my finite difference scheme but it did not occur to me to apply Maple's pdsolve() to it. It would have saved me some unnecessary work had I thought of it.  Thanks for pointing this out.

@wanderson The Fourier condition says that the heat flux is the same on both sides of x=L.  In fact, it's the purpose of the PDE to enforce that the heat flux is the same on both sides of any location; there is nothing special about x=L.  When we solve a single PDE on the entire domain 0 < x < 2*L, the Fourier condition does not enter the calculations because the PDE itself takes care of that.

But if you split the domain into two parts as you proposed to do in your initial post, the PDE on one side needs to communicate with the PDE on the other side.  That's where the Fourier condition comes in. Since I don't split the domain, that issue does not arise.