@sursumCorda 

With the linalg package the eigenvectors are returned this way

linalg:-eigenvectors(m);
                   [-1, 1, {Vector[row](2, {(1) = -1, (2) = 1})}], [7, 1, {Vector[row](2, {(1) = 1, (2) = 1})}]

But for the LinearAlgebra package you get a vector of eigenvalues and the matrix of eigenvectors.
If you set _EnvLinalg95 to true, use the option output='list' in  LinearAlgebra:-Eigenvectors (which then returns the eigenvectors under the form linalg-eigenvectors does):

_EnvLinalg95 := true:
LinearAlgebra:-Eigenvectors(m, output='list');
[[7, 1, {Vector(2, {(1) = 1, (2) = 1})}], [-1, 1, {Vector(2, {(1) = -1, (2) = 1})}]]

You have 3 indeterminates (x, y, z) and 2 equations (notice there is no f nor g as you say in the text).
Could you be more precise?

@Carl Love 

I tried hard, not enough obviously, to prove the equality by hand.

For the record the left hand side of rel(n) is a consequence of another relation, one can get on several papers, which is said to be "true for all integers larger than 3".
I observed it was also true for n=2, for any negative integer... and for any real number but 1 (so what is it to write the previous sentence?).
This puzzled me and then started simplifying rel(n), rewritting some terms differently and so on, thinking that assume/assuming would be of no help.
Thanks again for this Martin Gardner's Ahah intuition.

@acer @vv

Thanks to both us for your quick reply.

@Carl Love 

Sorry, I missed it.
But its notations remain confusing, not to say meaningless.

@Rouben Rostamian  

ContoursWithLabels is Kitonum's code see_here

@C_R 

Thanks for your return.

Go here help(dsolve[numeric]), search for 'known', try to understand hot to proceed from examples dsol8 and dsol8b.

@Muhammad Usman 

Does this suit you?

plots[display](
  PP
  , size = [500, 400]
  , axis[1]=[tickmarks=[seq(k=k/10., k=0..20, 5)]]
  , axis[2]=[tickmarks=[seq(k=k/10., k=0..20, 5)]]
);

But given L=1 and T=2, it seems to me that this

plots[display](
  PP
  , size = [500, 400]
  , axis[1]=[tickmarks=[seq(k=evalf(k*`&Delta:x`), k=0..20, 5)]]
  , axis[2]=[tickmarks=[seq(k=evalf(k*`&Delta:t`), k=0..20, 5)]]
);

would be more correct.

What is the structure of the data you want to send to Tecplot?
I guess it is a collection of contours, each contour being a collection of segments, that is a collection of 2x2 matrices.
The question is what structure Tecplot can manage? Once defined, the corresponding structure can be build in Maple.
But the "precision" of a contour of a given length is defined by the number of segments which approximate this contour. So the drawing in Tecplot of a contour made of, for instance, 100 segments constructed in Maple, will not have an higher resolution than thesame contour directly drawn in Maple.

As you defined

N := 20: 
Mx := 20: 

and plotted

plots:-contourplot(
        interfunc, 
        0 .. Mx, 0 .. N
        ...
)

why do you say "Here the t-axis is from 0 to 20 but its actual value is from 0 to 2 and the x-axis is from 0 to 20 but its actual value is from 0 to 2"?

If you want them to range from 0 to 2 just write

plots:-contourplot(
        interfunc, 
        0 .. 2, 0 .. 2
        ...
)

instead or set

N := 2: 
Mx := 2: 

It's up to you to lnow what you want to do.

About your claim "How I can, I extract data in the form of a dat file to plot in some professional software? " (which, I'm afraid of, will be likely understood as "Maple is not a professional software") you must be more detailed.
If you have professional in mind, I think that best thing to do is to give him interfunc plus some indications and let him manage to provide a "real professional layout".

@RaySierra 

And NO, I have no idea where your error came from.
What I can say is that your file, when opened with Maple 2015 on an iMac, generates a "Connection to the kernel lost" infamous message (which means that the only thing to do is to quit Maple).

@RaySierra 

There are probably hidden characters somewhere in your file that I don't have time to seek out (they could come from some copy/paste operation)
Just open the attached file and check if it works correctly with your Maple version
Test_cooling_task_model_rewritten.mw

Please keep me informed and let me know if you'd like an example of a Monte-Carlo simulation distributed over several cores (see my comment in my first reply).

Download Maple_PG2023_pp583-584.mw

About your Monte-Carlo simulatoin problem:
I use to do intensive MC simulations int he following specific context: some inputs X of a simulation code C are random variable and one is interested in assessing some statistics of the random output(s?) Y of C.

This problem is usually tackled by sampling X, running C over this sample, finally analyzing the induced Y sample. X sampling can be done through some strategy, in particular Monte-Carlo sampling (MCS).
MCS is largely used for a MC sample augmented by another MC sample is still a MC sample (specific properties MCS are conserved), opening the route to distributed computations, that is sampling+running on different cores/processors/machines.

This process is perfectly scalable meaning the sampling+running operation executed on P cores is exactly P times faster than executed on a single one (at least in theory for the OS of the machine host limits generally reduces the performance).
The time requires to gather the P partial simulations is general negligible

To do this I use to use the Grid package instead of the Threads one.

 @petit loup    @acer 

Sample(Normal(0, 1), [n, n])

@Jamie128 

So do you say that z = -2*M^2*sin(theta)^2 - 2*k*r^2?
What you need, I repeat myself, is a function r(theta, z) plot3d(r(theta, z), theta = a .. b, z = c .. d, coords = cylindrical)
and this expression of z doesn't make r uniquely defined:

restart
z = -2*M^2*sin(theta)^2 - 2*k*r^2;
r =~ [solve(%, r)]