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.

@Carl Love Your observation regarding the operation count of A^(-1) versus that of a triadiagonal solver is true.  In this case, however, the system A.U=B is solved m times within a for-loop with the same A, while B varies. In my mind I justified doing the inversion once and then applying the inverse m times as being not worse than solving the system m times, but I may be wrong on that.

@Chrono1 It is possible to generalize the technique to all other regular polyhedra. I will write up and post the general idea when I get some free time. For now, here is a demo of a geodesic on a cube.

@vv This is the best I could come up with.  Corresponds to the choice of parameters m=2, n=3, a=0.7.

@mmcdara Regarding "May we consider this exchange is over?":  Yes, sure.  I have nothing more to add.

@mmcdara I certainly did not aim to cast aspersions on your commendable attempts to explain, and if I have offended you in way, I do sincerely apologize.  My main point is that the solutions that you have exhibited so far are incorrect, and that's an objective statement, not a value judgement.  Here is what Maple produces for the solution to the null boundary condition case.  Note the clear depiction of the initial parabola and the zero values on the boundary. Compare that with yours.

restart;
pde := diff(u(x,t),t) = diff(u(x,t),x,x);
bc := u(0,t)=0, u(2,t)=0;
ic := u(x,0) = -x^2+1;
dsol := pdsolve(pde, {bc, ic}, numeric, spacestep=0.01);
dsol:-plot3d(x=0..2, t=0..0.5);

@mmcdara Your statement that the updator matrix is "always a quasi tridiagonal matrix" is not correct. If you formulate the problem properly, then the updator matrix will be exactly tridiagonal. This is done in many numerical analysis books.

In the correct formulation, the tridagonal updator matrix acts on the interior nodes only.  Boundary nodes (where the solution is known) are pushed to the system's right-hand side.

@mmcdara You are making this more complicated than it deserves, and along the way you are introducing errors.

To see what is wrong, set the boundary conditions and the source term to zero in your worksheet, as in:,
mu[1] := t -> 0;
mu[2] := t -> 0;
g := (t,x) -> 0;

The correct solution then should begin with the prescribed intitial condition and go to zero exponentially as time increases.  Looking at the graph of your solution, however, we see that's not the case.  Neither the initial condition nor the left boundary condition are satisfied, and the solution does not go to zero. 

@mmcdara I am afraid that there is something wrong with your code although I haven't gone over its details.

The code is attempting to solve the heat equation with a time-explicit finite difference method.  But that method will produce the correct result only if tau/h^2 < 1/2 in your notation. [That's called the Courant-Friedrichs-Lewy (CFL) condition.]  But in your case tau/h^2 = 15/2 which is far beyond the method's range of validity.

To see that something is indeed wrong, set the boundary conditions and the heat source to zero, that is,
mu[1] := t -> 0;
mu[2] := t -> 0;
g := (t,x)->0;

Then look at the solution Y at the end of your worksheet.  We see that Y correctly picks up the problem's initial condition but then it immediately drops to zero and remains zero at all future times.  It shouldn't be.  The correct solution of that problem will decay to zero exponentially.