@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);

@Pepini I see. Sorry, I misunderstood your wording.

What's your rationale for believing that an exact, explicit solution might be "seen" from the equations? I'm not seeing it (though my vision is poor).

I don't quite understand what's wrong with a numeric solution. Here a some ways:

restart;
kernelopts(version);
   Maple 2020.2, X86 64 LINUX, Nov 11 2020, Build ID 1502365

f := t -> exp(t) + exp(-t^2 - 2*t):
g := t -> exp(-t^2) + exp(-t^2 - 2*t):

S1 := fsolve( {f(a)=f(a+ba), g(a)=g(a+ba)}, {a,ba}, {a=-2..1, ba = -2..1} );
          S1 := {a = 0.1089038833, ba = -1.784296181}

S1 := fsolve( {f(a)=f(a+ba), g(a)=g(a+ba)}, {a, ba=-2}, {a=-2..1} );
          S1 := {a = 0.1089038833, ba = -1.784296181}

S1 := fsolve( {f(a)=f(b), g(a)=g(b)}, {a, b=-2}, {a=-2..1} );
           S1 := {a = 0.1089038833, b = -1.675392298}

plots:-display(
   plot([f(t), g(t), t = -2 .. 1], color=red),
   plots:-pointplot([seq(eval([f(a),g(a)], pt), pt=[S1])],
                    symbol=solidcircle, symbolsize=12, color=blue),
   size=[300, 300]);

@mmcdara I've already responded to that query about animations having a common range/view in the other thread.

It's an ancient weakness in the interfaces's animation plot drivers -- various so-called global plotting options are rendered common to all frames, for these kind of traditional animation. (Explore offers an alternative, frame-by-frame behaviour when animating.)

If I recall correctly, other plotting options with a similar drawback include at least lightmodel, tickmarks, gridlines, orientation, size, scaling, labels, axes. I'll note that for orientation a workaround is to use viewpoint (despite some inaccurate statements in the Notes section of it documentation).

@mmcdara Since you asked, the GUI renders a "traditional" animation by examining all frames and utilizing the largest range (view) necessary to include the features from all (together). And so all frames are rendered using a common range/view. By "traditional" I mean animations created with plots:-animate or plots:-display(...,insequence).

In contrast, the Explore command creates (and the GUI renders) the frames 1-by-1 on demand, and the rendering range/view of each frame is determined individually. Sometimes that is a bonus, as a way to work around the forced behavior of the traditional animation. Sometimes people want the fixed, common range/view when using Explore, in which case it can be supplied to the plotting command. (See my previous Reply's attachment, in which there is a commented-out call to Explore and I added a forcing common view option to the matrixplot call. You could remove that view option, to see the effect of each frame's tallest column determining that frame's vertical view.)