@Ronan To generalize my (GCD of numerators)/(LCM of denominators) technique, I think some restrictions need to be put on the field of the entries. Something like the field of quotients of a Euclidean domain. An example would be the field of rational functions with rational coefficients. Are your projective-geometry examples of that form?

@Aliocha Yes, degree 64 is not a practical example. I used it to emphasize the performance gain. AS I SAID AND YOU'VE IGNORED, I redid the example at degree 12 and I still got a 4X improvement. The number of test-data points used is irrelevant to any timing comparison, except that it needs to be sufficiently large so that neither timing is 0. Since the resolution of Maple's time is 1/64^th of a second on most processors (a huge value for modern processors), a large number of test points are needed to get over that hurdle.

I've occasionally seen processors for which Maple's time resolution was 1/256^th of a second, but I've never seen smaller than that.

@sursumCorda I think that the most likely explanation is that in those earlier versions, lprint and showstat didn't make use of screenwidth. The fact that the long number lprinted in your legacy version has no \ line continuation characters is good evidence that screenwidth is not being considered.

As an aside, long lines of code printed without hanging indents look hideous to me.

@Aliocha You are comparing two measurements that have been rounded down to 0; no valid comparison can be made.

You mention a "difference" in timings. Differences are irrelevant. It's the ratio of the timings that matters.

Please post a complete example.

My example shows a 9X performance boost for dense degree-64 polynomials. That's significant and substantial. You can lower the degree to 12, or whatever, and redo it. The performance gain may not be as much, but there will be some gain. For degree 12, I get a 4X improvement with Horner's.

@mmcdara Your P and H are identical because x is meaningless in this context. In other words, your converted polynomial is still in its original form, not Horner form.

@mmcdara Please post an error-free worksheet showing your timing test, and I will comment on it.

@Scot Gould Just eyeballing the OP's screenshot, I see that f1 and f3 are quadratic and the other four equations are linear wrt the 6 solution variables. That seems fairly simple, and it seems unlikely that Maple would choke on it. Because there are two quadratics, there are likely multiple solutions. So, you should change your eval to

eval(s, [T][1])

@RezaZanjirani Try this to get the side-by-side plots:

plots:-display(
    plots:-densityplot~(
        `<|>`(DiffrP||(1..3)), 0..1, 0..1, title=~ `<|>`("diffrP"||(1..3)),
        labels= [alpha, delta], labeldirections= ["horizontal", "vertical"],
        colorstyle= SHADING, restricttoranges
    )
);

There is some subtle and undocumented difference between 2021 and 2023 that caused the error. It's too minor to fuss about.

Regarding what you call "the traditional way": I don't see how plotting where certain functions are positive or negative can tell you which is the max.

Regarding multiple conditions: It's easy: Assign numbers, 1, 2, 3, etc., to the various combinations of conditions that you want. Then plot that number with densityplot.

@C_R I don't see any difference in the two strings you just posted.

I don't see any conjugation bar over any occurence of n or m in your worksheet. Please point out exactly where you see one.

The || opertaor allows even shorter syntax for sequencing:

`&#`||(945..969)||`;`;

@mmcdara mo, mi, and mn are needed when the name is part of a larger MathML structure, even if no modifiers are wanted.

How do you want to handle something like sin(ln(y))?

@MaPal93 I asked because if the solutions are equivalent, then it should be possible to verify that the first two are independent of gamma. One way is to take the derivative with respect to gamma and substitute a variety numeric values for the other variables, and gamma=1. Then use evalf, and you should get values close to 0.

@nm You could just use 

inttrans:-laplace~([ode_sys], t, s)