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

@Simon45 you write "I need to get a tridiagonal matrix somehow". The code that you have shown is very far from being there.  You need to fundamentally reformulate the algorithm. As vv has suggested, you should read your textbook to see what needs to be done.  Try to understand the algorithm. Don't think about Maple right now.  Maple will not help you to get there.

While you are reading your textbook, pay attention to what are called explicit method and implicit method for finite differences.  (Some books may call these the forward and backward methods.)

Have the following comments in mind as you read your textbook.

• The algorithm that you have shown implements an explicit (or forward) finite difference algorithm.  That method is only conditionally stable.  To obtain a correct solution out of it you need to make sure that the CFL stability condition is satisfied. (Read about the CFL condition in your textbook.)  In the code that you have presented the CFL condition is not satisfied, therefore even if you follow the suggestions made so far for fixing it, you will have a code that runs but produces junk!
• Once you have taken care of satisfying the CFL condition, you will find out that the solution does not deal with tridiagonal matrices at all.  That is, you are looking at the wrong algorithm if your aim is to arrive at a tridiagonal system.
• The tridiagonal system enters the picture when you do the time discretization through the implicit (or backward) finite differences. Study that method to see how the tridiagonal system comes into play.  By the way, you will learn that the implicit method is unconditionally stable (check your textbook) and therefore you need not worry about the CFL condition in that case.

@mmcdara Any function defined on an interval [0,T] may be extended periodically to the entire real line. You yourself noted this in your comment. So what's the problem with viewing the function f(t)=t^2 as periodic?

Note to the original poster: It is not clear whether you are asking for help with the mathematics of the problem, or how to solve it in Maple.

If you know how to formulate the problem mathematically, then it will be good to show that in your post so that people know that you understand what you are talking about.  Then many will be happy to show how to do it in Maple.

But if you don't know the underlying math, this is not the right place to learn it. Go back to your textbook or ask the teacher for help.

@Eric-pascal To get a sorted random vector, replace the line
X := LinearAlgebra:-RandomVector[row](n, generator=rand(100));
with
X := sort(LinearAlgebra:-RandomVector[row](n, generator=rand(100)));