@Paras31 My previous worksheet may have actually been a document in disguise (because although I removed all your 2D Input and replaced it with 1D input, it began its life as your uploaded document), and that may be why (totally guessing here) the plots didn't display. But the attachment below is a pure typed-from-scratch1D input worksheet.

Please download-and-open this worksheet directly in your Maple. Then execute it in its entirety by clicking on !!! on the toolbar. Then tell me if you get the same results as displayed below.

As always, plots will be crisper in your worksheet than the way they get rendered on MaplePrimes (which is why I usually copy-and-paste plots to MaplePrimes rather than post worksheets). In particular, the 3D plot in your "Task 2" will show fine detail rather than being a "blob".

Note that this worksheet computes and plots about 130,000 points in 6 plots all in under half a second.

Download OdeplotRefine.mw

@Kitonum I understand all that you did except for Step 4. How do you justify subtracting the angle from Pi?

@minhthien2016 Yes, I'm considering half planes. If the apex S is on the z-axis, and the base vertices are the vertices of the unit square (as we've all set up the problem), then the angle between any two faces that meet at an edge is at most Pi/2.

@minhthien2016 According to that Wikipedia page, a dihedral angle of a polyhedron is always an interior angle. I stand by the correctness of my Answer: Pi/3 and arccos(1/sqrt(3)).

@Paras31 Did you execute the code by downloading the worksheet from my Answer, or by copy-and-pasting the code from my Answer? I'm wondering if the problem is a 2D Input / Document Mode anomaly, things that I have very little experience with (or patience for).

@Andiguys The \n is the newline control character, I guess it doesn't work in 2D Input, which I never use. It's not a year-version issue, because I've been using it in Maple for over 20 years. It's purpose was to give extra space between the legend box and the w. You can just remove it: Change `w\n` to 'w'. Note that I also changed the quotes from backward to forward.

@emendes I'm not sure if you're asking a followup question, or you were about to and then figured it out. Does my selector (the one with And) work for you?

@Andiguys You could replace `πm` and `Πm` with pi__m and PI__m.

Here's how to do the plot (of course, there are a vast number of stylistic possibilities for presenting the same data):

W:= <seq(30..95)>:
TRC__min:= <seq(s[w][1], w= W)>/10^10:
Tau1__min:= <seq(eval(tau1, s[w][2]), w= W)>:
Tau2__min:= <seq(eval(tau2, s[w][2]), w= W)>:
plot(
    [<W | Tau1__min>, <W | Tau2__min>, <W | TRC__min>], 
    legend= [tau1__min, tau2__min, 'TRC__min'*10^``(-10)],
    labels= [`w\n`, ``], axes= boxed, thickness= 4
);

This plot uses the data from your original worksheet in this thread. If you've subsequently made changes that would change the numeric value of w, tau1, tau2, or TRC, then those new numbers aren't shown in this plot.

@emendes To select procedures whose names begin with d_, change procedure to 

And(procedure, suffixed('d_'))

in the select command.

@JAMET Are you incapable of making even the slightest generalization from my "advice"? The code in your Question contains 3 occurences of \\n. All three need to be changed to \n.

@nm The help page ?assuming says that (I'm paraphrasing) assuming a property in general without attaching it to specific subexpressions applies the property to all names in the expression. The problem is that y(x) is not considered a name, (but it is considered a function). This works:

simplify(residual) assuming positive, y(x)::positive

@JAMET This is one of your printf commands:

printf(`X=%a,   Y=%a   --->  %a\\n`, a, b, buff);
#This is the problem           ^^

Change \\ to \.

I am just reiterating @nm 's Answer (although I wouldn't say "\n is special").

In Maple 2023 on Windows 11, I've see A very frequently, B occasionally, and I don't recall ever seeing C.

@Andiguys I see in your wortksheet this procedure definition:

TRC:= (tau2, tau1)-> subs(DATA, `Πm`);

This is incorrect usage (although not necessarily an error) because the variables on the left side of -> do not explicitly appear on the right side, yet they appear indirectly via `Πm`. The indirect references will not be associated with the procedure parameters tau2 and tau1; instead, they will remain global variables. To create a procedure parameterized via indirect reference (this will also work for direct reference, but it's not necessary), use

TRC:= unapply(subs(DATA, `Πm`), [tau2, tau1]);

In Maple 2023+, MakeFunction is a synonym for unapply, because it makes more sense in common English.

@C_R I suspect that MapleSim doesn't have nontrivial arithmetic of matrices, polynomials, and even matrices of polynomials over finite fields, which is one of the uses of Maple's mod. By "nontrivial", I mean that the arithmetic includes determinants, inverses, and solutions of linear systems.

mod has two modes: modp and mods. The default is modp, which you're calling a "digital clock"; mods is the "analog clock". To change the mode, do `mod`:= mods before using mod in its usual form. `mod` is an environment variable, which means that if it's set in a procedure, it reverts after the procedure exits.

I don't understand why you object to using frem in the way you suggested.