@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.

@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.