@Ziaahmed812 For what it's worth, the following squeezes it a bit more (33% off the last timing above, for that size.):

restart;
Eratosthenes2a := proc(N::posint) 
local L, n, k; 
description "Calculate all primes less than or equal to N"; 
L := Array(2 .. N, i->true); 
for n from 2 to trunc(sqrt(N)) do 
  if L[n] = true then 
    for k from n^2 to N by n do 
      L[k] := false;
    end do;
  end if; 
end do; 
seq(`if`(L[n]=true,n,NULL), n=2..N);
end proc:

CodeTools:-Usage(Eratosthenes2a(2^20)):
memory used=132.76MiB, alloc change=20.00MiB,
cpu time=2.25s, real time=2.26s, gc time=988.78ms

nops([%]);

          82025

But I haven't looked at other accelerations.

@mmcdara I believe you would find that delaying evaluation (with uneval quotes, as you modified it) will not work in Malle 2021.1 for this example involving units. I had discussed this in my own Answer.

@Anthrazit I prefer to have a procedure do multiple-returns (ie, return an expression sequence), along with multiple-assignment. Eg,

  (A,B,C) := F(...);

with F here returning an expression sequence of three things.

That's for assignment/side-effects. For merely accessing globals I prefer using the colon-minus syntax for make references to a global variable, eg :-A to get at global A.

@GunnerMunk If I place the mouse-pointer on the output of the context-menu result (here, the result 9*Unit(m^2) from the inserted, inlined simplify call) and then make the keystrokes Ctrl= (`Ctrl` and `=` together) then that result of simplifying a*b gets updated with the new values for a and b.

In my Maple 2020.2 the consequence of pressing only the Enter key when the focus has progressed to the a*b 2D Input line is that I get a new, additional 2D Output (with updated values) , but the previous inlined context-menu result remains and is bumped down, as some kind of undesirable output artifact. That doesn't seem desirable, to me.

I don't use Document mode because I find the execution flow awkward, especially wrt to automatic advancing of the cursor when trying to use Enter or Tab, or lack thereof. I'm not telling you it's all ok. I'm simply informing you of ways that I was able to update the prior context-menu inlined output.

Your problems seem quite unspecific to units or the simplify command itself.

Your issues appears actually to be related to the execution flow that attains when on uses context-panel actions on 2D Input in a Document.

I do not know the "best" way to re-execute an inserted context-panel action (made on 2D Input, in a Document). But if I hit the menubar's !!! (triple exclamation button) then I can get all lines to re-execute sucessfully.

You did not attach your Document, which is not so helpful. I have guessed that it is like this:
  CM_Document.mw

note: I personally find Document mode's execution flow model to be unreliable and confusing, and only use Worksheet and 1D plaintext mode for input.

@Ronan  I see. The form you just gave didn't appear in your Question, but one way that it can be obtained is by using just a small part of the earlier code. Eg,

P := 14*c^4 + 84*c^3*d + 180*c^2*d^2 + 165*c*d^3 + 55*d^4 + 5*c^3
     + 21*c^2*d + 28*c*d^2 + 12*d^3 + 2*c^2 + 5*c*d + 3*d^2 + c + d + 1:

select(t->degree(t,[c,d])=3, P);

          3       2           2       3
       5 c  + 21 c  d + 28 c d  + 12 d 

If you need to guard against the case that polynomial P might consists of just a single term, then here is a reusable procedure that allows you to easily specify the degree and the variables,

F := (p,vars,deg)->`if`(p::`+`,select(t->degree(t,vars)=deg,p),
                        `if`(degree(P,vars)=deg,p,0)):

F(P, [c,d], 3);

           3       2           2       3
        5 c  + 21 c  d + 28 c d  + 12 d 

F(P, [c,d], 2);

             2              2
          2 c  + 5 c d + 3 d 

ps. The older link you gave is to a response from Carl Love, not me.

[edit] I made a correction to the code above.

@mmcdara I qualified my Answer with, "If I understand..."

I made no claim to have understood what the OP meant by his use of the word "extract".

The code I attached is due to an interpretation in which the terms are extracted -- eg. extracted separately into a list, rather than keeping the result as a sum of select terms with coefficients replaced/removed. (I believe that is consistent with how I interpreteted the word in an older Question -- if you're making a veiled allusion.)

I don't claim that I understand what is actually wanted here. The phrasing seems unclear.

It might even be that the OP wants to obtain programmatically the result,
    c^3+c^2d+c d^2+d^3
for the given example. And only a minor modification would be needed in the reusable/extensible code I attached, to handle that. As a shorter one-off,

P := 14*c^4 + 84*c^3*d + 180*c^2*d^2 + 165*c*d^3 + 55*d^4 + 5*c^3
     + 21*c^2*d + 28*c*d^2 + 12*d^3 + 2*c^2 + 5*c*d + 3*d^2 + c + d + 1:

`+`(map(t->t/coeffs(t,[c,d]),
        select(t->degree(t,[c,d])=3,`if`(P::`+`,[op(P)],[P])),[c,d])[]);

            3    2        2    3
           c  + c  d + c d  + d 

The use of operator form call `if`(P::`+`,...) above is to guard against the more general situation where polynomial P could be a single term instead of a sum of terms. The given example in the Question could also be handled with the following (if indeed this is what's wanted, and I make no claim that it is...),

map(t->t/coeffs(t,[c,d]),
    select(t->degree(t,[c,d])=3,P),[c,d]);

            3    2        2    3
           c  + c  d + c d  + d 

@JAMET If you haven't loaded any of the plots package then that ought to be,

  plots:-display(doPlot(a), doPlot(b), doPlot(rhs~(convert(c,list))));

and not merely (as you showed it),

  display(doPlot(a), doPlot(b), doPlot(rhs~(convert(c,list))));

which would return an unevaluated call to the (unassigned) name display, which would be an over-large output for printing.

@aksupriatna The CodeTools:-Usage command is used in the above code merely as a way of measuring how much time and memory the computation takes.

If your old Maple 13 version lacks that command then you can simply remove it from the above code. (The calls to the CodeTools:-Usage command merely wrap around the main computations.)

@srikantha087 If you want contours then why are you trying to use surfdata?

The surfdata command with its dimension=2 option does not allow for style=contour or style=surfacecontour. I have previously submitted requests/SCRs for that.

You didn't tag your Question as Maple 2019 when you asked it. (It's not so helpful to us, to leave out such key information. I've tagged it now, myself.)

In older Maple 2019 the contourplot command doesn't support a legend for the contours. It also doesn't support colorscheme.

[edit]
The listcontplot command also doesn't support legends or colorscheme, but after using Interpolate:-SplineInterpolation on the Matrix of data then we can obtain a bivariate expression in lieu of discrete data.

In this old answer I implemented a way to use colorscheme with a 2D contour-plot. And in this old post I implemented a way to contruct a colored legend for a 2D contout-plot. I've never merged together those two pieces of added functionality.

So,... if you want a legend of colored items (akin to a finite number of contours) along with a 2D density-plot then see my Answer below.

But if you really want a 2D contour-plot (with discrete shading) that allows for your multi-color colorscheme as well as legend items then please indicate that clearly.

@mmcdara Yes, it really is wrong, because Tom used the claim (about surfdata's producing only 3D plots) to support the idea that -- in consequence -- they do not allow for legends.

The point is that only 2D plots can have legends, in the technical/built-in sense. 3D plots can have titles and captions, and those may be used as a kind of workaround for a legend, but without functionality to specify location. 3D plots don't allow for the legend plotting substructure.

@tomleslie Your claim about the result from a call to surfdata being a 3D plot is not true. The plots:-surfdata command can produce either a 3D surface plot or a 2D density-style plot.

And the calling sequence example that the OP provided in his Question shows the calling sequence for the 2D result.

Please provide your Matrix M, or full code to reproduce. That can be done using an uploaded and attached worksheet (green up-arrow in the Mapleprimes editor).

Also, you should state explicitly what you want.

One possibility is a "colorbar", ie. a thin bar with smoothly shaded gradient, normally displayed above or beside the plot. That's possible, but tricky if you want the whole assembly arbitrarily sized or if you want to be able to export it as a single image.

An easier possibility is a discrete number of shaded items (lines or thicker bars, corresponding to discrete heights or contours) displayed using the usual legend mechanism. If that's what you want then how many such items do you want? Do you want a fixed number, just for your fixed example with its particular kind of colorscheme? Or do you want re-usable and flexible code that can handle arbitrary numbers of such legend items for arbitrary coloring schemes?

The kernel (Maple engine) builtin arithmetic operators `+`, `*`, etc, are not units-aware. As far as they are concerned units are just some extra function calls which some expressions just happen to have as multiplicands.

Hence resolving and combining units can be done by Maple Library level commands. That's shown in my Answer's latter group of examples.

The Units:-Simple and Unit:-Standard packages can be "loaded", which rebinds various command names (including `+`, `*`, `/`, etc) with Library-level versions which are units-aware. These versions can resolve and combine units automatically. These were my Answer's first groups of examples. In my Answer I included URL references to the Help pages of those packages, so that you could read more about them.

@Felipe_123 

restart:

f:= proc(p,c) local j;
      Matrix([seq([0$j,
                   PolynomialTools:-CoefficientList(p,x)[],
                   0$(c-degree(p)-1-j)], j=0..c-degree(p)-1)]);
end proc:

f(x^4+3*x^3+x^2+2*x+1, 10);

f(x^5+11*x^2-3*x, 10);