@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)));

tizozadoxo, it is understandable that you wish to pass your course, but having other people do the work and then handing that work to your teacher as if it were your own, is not the best approach.

At least, attempt to solve the problems yourself. If you run into difficulties, show your work here and ask for help on specific issues that you are having. But just posting the statement of a homework problem and asking for someone to solve it for you is not quite the right thing.

@tizozadoxo Here is a stripped down version. I can's say that it's "easy" because that depends on what you already know in mathematics and Maple.

restart;
frames := seq(plot(1/(1+e*sin(t)), t=0..2*Pi, coords=polar), e=0..2, 0.1):
plots:-display([frames], insequence, scaling=constrained, view=[-3..3, -5..2]);

I have written up a detailed analysis of the rolling of an ellipse along a straight line.  I will put it under Posts in MaplePrimes, so you may want to check it out.

I haven't looked at the question that you have asked, that is, rolling of an ellipse on anything other than a straight line, but I expect it to be significantly more complex.

You should explain:

• Which boundary layer equations are you interested in?
• What is the "hmp method"?
• Why do you want to solve with the hmp method rather than by some other method?