I would not use remove(has, ..., unit) since that would act quite differently on expressions other than a product. (Your expression assigned to a is of type `*`, but it's unclear that is all you want to be able to handle robustly.)

I would suggest these preliminary steps be part of the process:
1) combine addends into common units
2) optionally convert to a preferred unit (or accept the base unit for the current system, default=SI)
3) use convert(...,unit_free) to strip off the remaning units, as suitable

For example, quoting unprotected global versions of option names for robustness,

Download remove_units.mw

Maple was originally designed for efficient symbolic computation, as opposed to pretty mathematical typesetting.

A key aspect of that original design was the storing in a session (by the Maple engine, a.k.a. "kernel") of only one representation of any sum or product of terms (commuting under usual addition and multiplication). Any other input of an equivalent form of an expression -- with terms commuted -- gets uniquified to the stored form.

Hence, if the form
   -b + a
is first formed in a session then construction/entry of the alternate form,
  a - b
gets automatically recognized as the earlier, stored form.

However, for this kind of expression the terms can be reordered by the sort command. For example,

Download simpl_examp.mw

As far as mathematical typesetting goes, one could also display the expression (in marked up "2D" math), without the above issues interfering. Here I use so-called inert operators, so avoid those uniquification issues,

Download inert_add_ex.mw

If you have other, particular usage scenarios then you should describe them precisely, as a best bet on getting adequate assistance.

You asked about plotting the elements of a sequence, so I'll address that aspect.

You can use the plot command to show a discrete collection of x-y points.

The data can be formed in a list of lists, or a Matrix. (I'll do it as a list of lists, below.)

You can adjust other qualities, such as the symbol, the style, the color, etc. See the Help page plot,coptions for more details.

Download point_plot_ex.mw

The plot shows that the local maximum is the global maximum over that range t=0..2.

So you can use the Optimization:-Maximize command.

In my opinion this is more straightforward, flexible, and likely more accurate in general than either examining the discrete data points from a plot structure, or using an interpolated spline.

I use the output=listprocedure option for dsolve, as a convenient way to get at procedures for xn(t) and zn(t), similarly to what I showed in your previous Question.

Download evalv_max_ac.mw

There's nothing intrinsically wrong with returning NULL as an explicit statement. This is a common programming technique in Maple.

   return NULL;

One alternative is to have a statement that does nothing and happens to also return NULL, but that would be obfuscation for no benefit. The explicit return call is clear and understandable programming.

ps. Using assign brings the burden of having to remember that happens to return NULL. Since you didn't already think of that then using it would entail one more thing for you to remember. You already know what explicitly returning NULL means, and it's more likely that others reading your code would too.

Without a range supplied fsolve is trying to find roots from A=-infinity to A=infinity. That search is missing values very close to zero.

[edit] The situation is made more difficult due to the root lying close to a singularity, where it is steep, and with the expression nonreal on its other side.

Download A_ac.mw

Here are two ways:
1) using plots:-odeplot directly with the sol1 returned by dsolve.
2) using plots:-spacecurve with procedures returned in sol1 by dsolve.

The second way requires some use of eval to extract the various results (after applying procedures at a numeric t value). There are a couple of ways to do that. I show one way that utilizes the output=listprocedure option. It can also be done without that option, using eval slightly differently. The first way looks simpler to me.

Download unable_to_evaluate_ac.mw

@Syed Asadullah Shah

Your code doesn't initialize F (or T) as an Array, so your assignments to indexed references cause a table to be assigned to F (and T). But your code also shows F(k+2) , which is nonsense. What do you intend by that?

Your attempted second call to pade (inside the solve call) is messed up 2D Input, and is actually being parsed as pade*(...) . You could get rid of that multiplication.

The n that appears in the expressions returned by your procedure gnrpoly is the local variable n of that same procedure.

That escaped local name n is not the same as the global name n, which you are trying to use for substitution at the top-level. That is why your calls to eval are not making a substitution.

There are several ways to fix your issue.

One approach is to utilize unapply (w.r.t that local name n) inside gnrpoly, and have gnrpoly return a procedure. Another way is to pass in the variable name (eg. n, from the higher level) as another parameter of gnrpoly. Below I show the former approach, with unapply.

Download Q_10-10-22_gentrate_eqn_ac.mw

ps. You used the term "inert". This is not an issue of inertness.

I disagree with your claim that "foo() has z" used.

The only z that foo "has" (apart from its declaration as a local of foo) is as a parameter of an anonymous procedure defined and used within foo. And that is not the same z as the local z of foo.

The local z of foo is not used within foo.

An alternative is,

  plots:-textplot([2,2,InertForm:-Display(%floor(5.5),'inert'=false)])

You could use the inert form %log[b] . You might also use InertForm:-Display(..,inert=false) to get that rendered in black (instead of gray).

DisplayLogQCM_ac.mw

How about specifiying the types on the (procedural) parameters in the preliminary call to codegen[makeproc]?

Eg, (notice the parameter specifications of the procedure),

   codegen[makeproc](TheList, parameters = map(`::`, [proc_params], float[8]))

Those were specified programmatically, not manually.

Download C_code_generation_with_CodeGeneration_ac.mw

Your current problem is due to the fact that your collatz procedure is returning NULL.

That is, calls like collatz(13) have a return value of NULL. That happens because the last line of your collatz procedure merely prints max_i. The print command itself returns NULL.

Try it with the very last line of collatz being,

   return max_i;

after (or instead of),

   print(max_i);

Here's one simple way,

   m:=Matrix([[1,2,5,6,9],[3,4,5,13,14],[1,3,4,10,11]]):
   dataplot(m^%T, style=line, gridlines=true, size=[500,"golden"]);

The size option is optional. (That instance of "golden" refers to the aspect ratio; I could have used approximately 0.618 instead.)