@salim-barzani OK, lets take epsilon=1 then I agree with Eq 3.7. I needed W because you can't integrate wrt a function in Maple, i.e. int(..., w(xi__1)) doesn't work.

But as I showed there is no way to transform F to F2, so the integral that gives the arctan result is invalid, and so Eq 3.13 is incorrect, so it will not satisfy the ode. You can perfom the integral with e2 and e0 in (just change the Int to int). The problem is then that you get xi__1 = f(w) but here and in general f(w) has multiple w in, so you cannot easily invert it to get w as a function of xi__1.

The solution to this is to solve the ode in Eq 3.7. dsolve will do this directly, or you could take the square root and integrate. Since you get a general solution it will work for arbitrary e2 and e0. If you want to then find cases that are real for various choices of e2 and e0 you can make assumptions that make the things under the square roots positive.

@salim-barzani Almost there.

@salim-barzani Well if it were changed I might look at it. I have several times pointed out that it makes it harder to help you, and I have ignored several of your recent questions for this reason. In the present case using w instead of r will be especially difficult since the paper uses w later for a different purpose.

I notice that as usual you changed the notation in your worksheet to make it harder to follow the paper. Perhaps you wish to correct that.

@sand15 Sorry; I must have google searched it instead of clicking the link

@sand15 Your second link does not work (not a web address).

@Alfred_F It is a pity that the nice numbers aren't really proven because the maximization works in floating point numbers, so r = 0.749999999968849584 we guess is 3/4, and identify guesses at the exact values of the angles. Possibly by doing the differentiations involved in the gradients involved in the maximization and solving, you could do it all symbolically.

geom3d is a nice package and easily works here for the duals since the radius of the circumsphere and in-sphere are built-in. It has a duality command that I didn't make use of. One limitation, unlike its 2D equivalent geometry, is that the sizes have to be fixed at the time of creation, so you can't have a cube with sides a, where a is a name.

@dharr Yes it works with 3 Euler angles; see updated answer. I should have realized that right away. As usual you have to play around with the initial point.

@Alfred_F If I do the further tilt you describe, I can get 9/16.

cubeoctOptimization2.mw

I suspect I need 3 Euler rotations, so I'll try that and update if it works.

Can't be sure it is a global maximum. [Edit: updated to use Euler angles, which gives a larger maximum than my first attempt.]

Download cubeoctOptimization2.mw

Even though kernelopts(version) is incorrect, interface(version) correctly reports 2026.1.

@C_R Converting the lprint(%) to 1-D resolves the issue. But I was not able to generate it with a fresh 2026.1 worksheet (I didn't try document mode). But I have a recollection that recently 2D math with % in has not worked consistently. So it might not be to do wth lprint but with % in 2D. Unfortunately I can't get it to not work when I try a few examples. If I come up with one I'll post it here.

You added a second question with the uploaded file instead of editing your own question and uploading the file to it. So I did this for you, moving your file to where it should have been. I then deleted the now redundant question, with the reason "duplicate question" and not "spam", though perhaps the system shows you that.

I am happy to accept donations for my brain upgrade. On the other hand you may wish to delete this post.

Maple's standard diff can differentiate with respect to theta, but not to theta(s). However diff is modified to accept theta(s) when the Physics package is loaded. So with(Physics) solves the error and returns some output, but I don't know if it is what you wanted.

@Mth2026 That download page for the assistant implies that it is not included with Maple Flow but has to be downloaded separately.