@imparter You have shown the expected graphs of w(r) at z=1.  To get those graphs, set z=1 in your formula for h, calculate h, then use that h to plot the graphs as I shown in my earlier message.

By the way, in my code I have "M in [2, 5, 7 ,9]".   That's not correct.  The numbers in the square brackets are expected to be the values of M^2, not M.  To correct that, replace "M in [2, 5, 7 ,9]" with "M in sqrt~([2, 5, 7 ,9])".

@imparter I do realize that h is not a constant.  That's why I presented two samples with h=1 and h=0.8.  You may pick any h that you wish.  Or if you prefer, pick a value of z first, and then use your formula to calculate the corresponding value of h, and then use that h to plot the graphs.

I don't understand what you mean by "how to combine all those in one graphs for different range".  What function do you want to plot?

@imparter 

Domain is 0 to 1?  You must be kidding  :-(

Your unknowns w(r,z) and theta(r,z) are functions of two variables, r and z.  So your domain is a region in the two dimensional r,z plane.

To set up the PDEs properly, you need to

1. specify the geometric shape of that two-dimensional region;
2.  identify that region's boundary (precisely, not just hand-waving);
3.  specify the conditions that the unknwons w and theta take on that boundary.

@imparter You are attempting to solve a pair of PDEs in the unknowns w(r,z) and theta(r,z)..  You specify the boundary conditions as

conds:= w(h(z),0)=0, theta(h(z),0)=0, D[1](w)(r,0)=0, D[1](theta)(r,0)=0;

What is the domain of the functions w(r,z) and theta(r,z)?  You should specify the domain clearly and unambiguously before attempting to prescribe a meaningful set of boundary conditions.

1. The definition of h involves delta.  You need to assign a numerical value to delta.

2.  In the definition of pdes you have (DP+theta*x4*sin(alpha)).  Shouldn't that theta be theta(r,z)?

@one man and @dharr 

In this revised worksheet I have incorporated the improvements suggested by dharr which reduce the time taken by dsolve() from four minutes to 2 seconds(!) and also extra code to calculate the length of the geodesic curve.

Additionally:

1. I have changed the independent variable in the geodesic's parametrization from the previous u to v.  Now dsolve() works with no extra options.
2. As noted by one man, there is a second geodesic.  I have added code to calculate the second geodesic which actually is a local maximum of the arc length.  The lengths of the two geodesics agree with those reported by one man.

Here's what the two geodesics look like.  Their lengths are 3.7262 (black curve) and 3.8439.(magenta curve).

Download geodesic-on-rounded-cube2.mw

@goebeld Yes, that's also possible.  And if you are going to load the Student:-MultivariateCalculus package, you might as well calculate the norm through Norm(T).  The advantage of the "^+" operator is that it's always avaialble and does not require loading a package.

@goebeld the "^+" applied to a vector or matrix produces the object's transpose.  Therefore, if A and B are column vectors, then A^+ . B is their dot product.  In particular, sqrt(A^+ . A) is the length (that is the Euclidean norm) of the vector A. 

Alternatively, the dot product and norm can be calculated through dedicated functions in the LinearAlgebra package, as in
LinearAlgebra:-DotProduct(A,B, conjugate=false);
LinearAlgebra:-Norm(A, Euclidean, conjugate=false);
These are too verbose for my taste.  I prefer the "^+" notation.

Aside: There is also the older "^%T" operator for producing the transpose, as in A^%T.  but I have no reason for prefering that over the shorter "^+".

You will need to be more specific regarding "angle in the x-z plane".   That can be (a) angle relative to the x axis, (b) angle relative to the z axis, (c) angle relative to the xz plane.  It also can mean the angle relative to the x or z axes of the vector's projection onto the xz plane.  

There are possibilities beyond those five.  For instance, you may specify a vector through its spherical coordinates (azimuth angle, polar angle, and length).

To remove any ambiguity, it will help if you could provide a concrete example of how you would like to describe a vector by its length and angles.

1. The title says "3d" but the question says "2d".  Which is it?
2. What do you mean by "edge coordinates"?
3. What does it mean to "plot a coordinate"?

@Christian Wolinski That's a very clever combination of map and map2, and does the job well.  Thanks!

@mmcdara You wrote: "The case of spherical coordinates is more complex (read the reference I provided you in y previous reply)".  

In the reference that you have provided, integration in spherical coordinates is made unnecessarily too complicated due to the attempt to integrate with respect to radial direction first and then with respect to the angle.  The solution is simpler if we integrate with respect to the angle first.

@CaptNnFalcon Your English is fine, and being new to Maple is not unusual.  The starting point, hower, is to draw the diagram that I suggested.  Do it by hand, not in Maple.  You can't get started without drawing that diagram.

As this is a homework problem, you are expected to do the work rather than asking others to do it for you.  At least provide some indication of the effort that you have put in toward solving this problem.

Suggestion: If you haven't yet done so, it is essential to draw a good drawing of the intersection of the torus with the y=0 plane.  Mark on the diagram the coordinates of the intersection points with the axes.

@acer I see no point in attempting to optimize OP's "dummy code" (his words) as we don't know what his real application looks like.

But yes, if we want to optimize the dummy code, we call fun just once within the proc, as in F := fun(x,y), and then use F as many times as needed.