@zenterix How is that person executing the .mpl file? By running command line Maple? Or (as you suggest above) by reading it into a worksheet? The way I generate packages or help databases is by running a worksheet which gets other files as necessary from the same directory, and uses worksheetdir. (Getting them from a subdirectory is a simple modification of that.) Someone else can run that worksheet which is in the same directory as all the required files on their system, since I would give then the zipped contents of the directory.

If you read a .mpl file with read, then the code is running in the worksheet you read it into, not in the .mpl file.

@mmcdara Yes, I was developing some code like that before I realized I was reinventing the wheel. You can use generator=(()->1) so you always get 1, then density gives the probability. Actually not quite because the ones along the diagonal were removed. Since the native graph datastructure has to be made from the matrix, this may not be as efficient as directly generating the edges and vertices.

@zenterix That's why I'm confused. You are opening a worksheet to read the .mpl file. So in order to read the .mpl file at this point you know where it is. So what do you mean by getting the path "programatically from within the mpl file"?

Note that once the worksheet is saved and then you exit Maple and reopen, the worksheet is at currentdir(). So if I distribute some code in a series of files, I can distribute a worksheet in the same directory that reads test.mpl or does something else. Then when the user runs that worksheet, currentdir() (or worksheetdir) will be where that worksheet is and where test.mpl is.

@JohannesC My Maple version (2015) had some errors so I cut and pasted CONV_2 to a new worksheet. But now you have inserted the missing line, I am getting the same result as you. 

The problem is that rho has an internal form that is different from the way it looks (try lprint(rho) to see this). The workaround is to write in the expression for rho in the line that generates CONV_2 - then simplify works.
 

Download Tensor_Calc_simplified2.mw

@WA573 
 

Download CM_TW2.mw

@zenterix I added an example to my answer, using LibraryTools as @acer pointed out.

@AHSAN So you will need to arrange your loops so that the value of Q and lambda are known at the time you call dsolve.

An outer loop can change values of Q and lambda and the innner loop with the other parameters as you have now.

@AHSAN Suggest you fix it then.

You define the function theta(y) in terms of y. Then you use dsolve to find u(y), but you give boundary conditions for the known function theta(y), not the unkown function u(y).

@moh111 What do you know about the signs or ranges of k, re and the variables?

@tomleslie The OP specified the 8 unknowns they were interested in. My understanding is that the OP wants these as functions of k and re. Setting infolevel[solve]:=5 shows that solve uses methods for polynomial systems, but they are very slow. (If k and re are specified then fsolve easily gives a solution.)

You should not use both N(x,t) and N[0] together - use N__0 instead of N[0] (looks the same), and the same for omega, E etc.

You set a Matrix equal to scalar zero; I assume that was meant to be the determinant of that Matrix.

If I guess at what substitution you want, then I get the attached, which still has cosines and x and t in. You will need to specify exactly the steps that lead to the result you expect. My guess is that it is not a matter of substitution, but of solving a PDE and getting some eigenvalues for E, but the equations aren't familiar to me, so many more details are required.

d_relation.mw

@bstuan I was just too lazy to write out the 2023 case :-)

Edit: That is, (A.B)^n is not equal to A^n.B^n in general, though there are some cases where it is true e.g., commuting matrices, n=1, n=0.

@bstuan Note (A.B)^2 = A.B.A.B <> A.A.B.B

@Earl It's not clear to me why the loxodromic spiral would be the (2,3) case. There's only a single spiral, so it's more like (2,1), but that doesn't work out quantitatively. Here is separate code for the loxodromic case:

Download LoxodromicCircles.mw