Your procedure extractxishu contains the command,

    remove(has,i,{p,q})

But the procedure extractxishu has the name p as its parameter. Therefore the "p" inside that call to remove is a reference to the parameter p of the procedure, and not a reference to the name  p that appears inside your polynomial (ie. the polynomial that you pass in to the procedure). That's the fundamental mistake and reason why the "p" appearing in the polynomial is not being removed by the procedure call.

One way you could correct that is by making the "p" inside that call to remove be a reference to the global name p that appears in your polynomials. Eg,

extractxishu:=proc(p)
local i,xishulist,teamlist;
  teamlist:=convert(p,list);
  xishulist:=[];
  for i in teamlist do
    if type(i,`^`) then
      xishulist:=[xishulist[],1];
    else
      xishulist:=[xishulist[],remove(has,i,{:-p,:-q})];
    end if;
  end do:
end proc:

Another simple correction would be to change the name of the procedure's parameter (and the places it gets used, such as the conversion to list), to avoid the collision. That could allow your procedure to still pick up the p and q reference, scoped  from the higher level. I haven't shown that, but it's a simple edit.

Or, you could pass the set of names in as as an additional argument to the procedure. Eg,

extractxishu_alt:=proc(p, varset)
local i,xishulist,teamlist;
  teamlist:=convert(p,list);
  xishulist:=[];
  for i in teamlist do
    if type(i,`^`) then
      xishulist:=[xishulist[],1];
    else
      xishulist:=[xishulist[],remove(has,i,varset)];
    end if;
  end do:
end proc:

extractxishu_alt(eq2, {p,q});

I'm not convinced that your procedure would always produce the result you'd expect, for additional examples, partly because I don't understand all your description. But you've already indicated that you are primarily interested in knowing what went fundamentally wrong with your original, so I've tried to focus on that aspect.

The error comes from premature evaluation of the first argument passed to the plot command, ie.,

funktion(x)
Error, (in funktion) cannot determine if this expression is true or false: 2 < x*Units:-Unit(1/m)

That evaluation is part of Maple's usual evaluation model for procedure calls. The procedure funktion is not set up to deal with the symbolic name x as argument, and it is being called before variable x takes on any numeric plotting values.

Using Maple 2021.1, we can work around that by using the so-called operator form calling sequence of the plot command, eg.

restart;
with(Units[Simple]):
funktion := x -> if 2*Unit('m') < x then 1/x; else 2*Unit('m'); end if:
r1 := 3.0*Unit('m'):
r2 := 6*Unit('m'):

plot(funktion, r1 .. r2);

Normally we could also work around the issue by alternatively: 1) delaying evaluation (single right-ticks around the first argument to plot) or, 2) handling the premature evalution by having the procedure return unevaluated if its first argument if not as expected. Those alternatives are problematic here, because of the way that the plot command currently (Maple 2021.1) handles units. (I've run across that before, and consider that a bug, which I've previously reported. It ought to be able to handle those alternatives just as well as it does the operator calling sequence.)

The result of parse("l") is the global instance of that name, and does not match the local l of your procedure.

So, in your procedure, when loop index j is "l" then dummy1 gets assigned according to the global name l and not according to the local name l.

Here is a simpler procedure to illustrate, without the looping or embedded component references.

restart;

test := proc()
  global G;
  local L;
  G := 3.2;
  L := 5.7;
  print("G", eval(parse("G")), G, :-G);
  print("L", eval(parse("L")), L, :-L);
end proc:

test();

        "G", 3.2, 3.2, 3.2

        "L", L, 5.7, L

L := 99;
                               99
test();

        "G", 3.2, 3.2, 3.2

        "L", 99, 5.7, 99

In my opinion using declared globals in a procedure is generally a poor programming practice. Using parse to turn strings into names -- within a procedure -- is also generally poor a programming practice. I think that these should be avoided where feasible to do so. Even large and involved applications with many embedded components (even a dynamically generated collection of same) can often be programmed without resorting to such technique.

You are encoutering premature evaluation of,
   diff(coeff(5*x+x*y+4*x*y^2,y,k),x)
before k attains its numeric values.

One alternative is to use add instead of sum, where the former has special-evaluation rules that delays the evaluation until k gets the numeric values. I generally prefer that instead of using single right-quotes (uneval quotes) to delay the evaluation, if I'm just trying to add up a finite number of items and not doing actual symbolic summation.

In the example below I replaced xy^2 with x*y^2, for fun.

restart;
c:=proc(P) add(diff(coeff(P,y,k),x),k=1..2);end proc:

c(5*x+x*y+4*x*y^2);

            5

restart;
c:=proc(P) sum('diff(coeff(P,y,k),x)',k=1..2);end proc:

c(5*x+x*y+4*x*y^2);

            5

note. This is a common usage misconception. It is essentially the same underlying issue as in your Question from September 2019.

The output=XML makes the Tabulate command return a call to DocumentTools:-Layout:-Table.  That can be used in other Layout:-Table instances, attaining a nested effect.

Here is a (somewhat uninspired) nested example,

restart;
with(DocumentTools): with(DocumentTools:-Layout):
T1 := Vector(5,i->i):
T2 := Vector(3,i->plot(x^i,x=-1..1,size=[100,100])):
xml1 := Tabulate(T1,output=XML):
xml2 := Tabulate(T2,output=XML):
InsertContent(Worksheet(Table(exterior=none,interior=none,
                              Column()$2, alignment=center,
                              widthmode=pixels, width=200,
                              Row(Cell(xml1),Cell(xml2),
                                  alignment=center)))):

nested_Tabluate_ex.mw

Also, I have seen a few cases where the Tabulate command didn't offer some esoteric finer control that I wanted, and after using output=XML I was able to programmatically adjust the raw result and then embed it using a call like the above, ie.,
   InsertContent(Worksheet(... my_xml ...))
Other related use of output=XML is to be able to save the result to a new .mw file, or to embed the result in a new worksheet tab, eg.,

restart;
with(DocumentTools):
L := [seq(plot(x^i, x=-1..1, labels=["",""],
               size=[200,200]), i=1..3)]:
xml := Layout:-Worksheet(
         Tabulate(L, output=XML,
                  widthmode=pixels, width=600)):
:-Worksheet:-Display(ContentToString(xml));

If I understand your question this might get you started.

In the following attachment I create an appliable module (for minor efficiency, so that the inner proc is not created each time F gets called). But a simpler revision, which still allows reusability wrt the degree and variables, is,

F:=proc(p,varlist,deg)
  map(proc(t,vars) local u,v; u:=coeffs(t, vars, v); [u,v]; end proc,
  select(t->degree(t, varlist)=deg,`if`(p::`+`, [op(p)], [p])),
        varlist); end proc:

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:

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 ]]

cidea.mw

I supposed that you'd want the coefficients associated with the corresponding terms, hence the returned lists of pairs.

You could also get just the coefficients, though naturally then you wouldn't know which terms contained them. Eg, with P being the polynomial,

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

          [5, 21, 28, 12]

Your mistake lies in your passing an expression sequence to the CodeGeneration:-Matlab command. That makes that command treat the second argument as if it were an (invalid) option to that command.

The result from solve, for your example, is an expression sequence of two scalar expressions.

To remedy this, you could call Matlab(...) on them separately, or put them into a list before the call, etc.

By the way, it seems that the expressions could be simplified and made more compact.

See attached.

calc_theta_ac.mw

ps. You had your Question marked as "Maple 18" (a version from the year 2014). But your attachment was last saved in Maple 2018, which is a later version. I have adjusted the Questions marking accordingly.

If the name begins with an ampersand then the parser will accept its use as a binary infix function call (ie. allow its use as a binary infix operator).

For example,

restart;

`&A` := (a,b)->a^b:

x &A y;

               y
              x 

You don't have to assign anything to &A, to construct the call.

Also, you can also enter the call using prefix form, ie.  `&A`(x,y) , and by default that pretty-prints in the infix form.

The above doesn't allow you to enter your stated, literal example   x A y  as input, without the ampersand. However you can get some interesting effects for pretty-printed output and computation, by using the so-called extension mechanism of some low-level commands. Eg,

restart;
`print/&A` := proc(a,b) uses Typesetting;
                mrow(Typeset(a),
                     mi(" A ",':-fontweight'="bold"),
                     Typeset(b));
              end proc:
`expand/&A` := () -> A(args):
`value/&A` := () -> A(args):

A := (a,b) -> a^b:

foo := sin(x) &A sqrt(y);

expand(foo);
value(foo);

In order to alleviate burden on the GUI plot renderer the result of HeatMap is usually accomplished (internally) by use of a background image.

In Maple 2021.1 the background image is being incorrectly rendered to fill the whole inlined plotting window, instead of (as previously) filling only the area bounded by the axes. (My Maple 2021.0 for Linux does not seem to have the problem, so it might be new to point-release 2021.1.)

Two possible workarounds are (undocumented) forcing of the result as either a collection of polygons or a densityplot.

Those put more burden on the GUI, however (ie. they render more slowly, and can make the GUI response sluggish if the plot is swept-selected or right-click, and in some cases when the sheet is scrolled). Export to PNG format image seems reasonably sharp.

Statistics:-HeatMap(Matrix(128,(i,j)->modp(binomial(i,j),2)),
                    color=["White","Black"],method=polygons,
                    axis[1]=[location=high],axis[2]=[location=low],
                    size=[600,400]);

Statistics:-HeatMap(Matrix(128,(i,j)->1-modp(binomial(128-j+1,i),2)),
                    color=["white","Black"],
                    method=densityplot,style=surface,
                    axis[1]=[location=high,tickmarks=[seq(i=i,i=1..128)]],
                    axis[2]=[location=low,tickmarks=[seq(128-j+1=j,j=1..128)]],
                    axes=frame, size=[600,400]);

Here are some ideas, which might help you get started,

Shane_Kreller_M6_Maple_Assignment_ac.mw

Here is another way, using seq instead of a loop,

Given you initial list assigned to z,

LLq:=iquo(nops(z),26):
LLr:=irem(nops(z),26):
AZ:=parse~([$"A".."Z"]):
ans:=[seq(seq(`if`(i=1,AZ[j],cat(AZ[j],i-1))=z[26*(i-1)+j],j=1..26),i=1..LLq),
      seq(`if`(LLq=0,AZ[j],cat(AZ[j],LLq))=z[26*(LLq)+j],j=1..LLr)];

You might also create AZ with,

AZ:=convert~([$"A".."Z"],name):

Mapleprimes_Label_lists_ac.mw

Naturally, it's also possible to start the LHSs with A1..Z1, etc,

[seq(seq(cat(AZ[j], i) = z[26*(i-1)+j], j = 1 .. 26), i = 1 .. LLq),
 seq(cat(AZ[j], LLq+1) = z[26*LLq+j], j = 1 .. LLr)];

Here is one way to construct a 2D surfdata plot that has a legend.

The code is implemented as a procedure that accepts your original surfdata arguments, your simple colorscheme type (a simple list of colors), and also a few extra options to specify the number of legend items, legendstyle (location), etc.

It's a bit rough, but might serve your purpose. With a little effort I could make it more flexible.

Here is an example with data I made up (because you hadn't yet provided your data.) [edit] I correct an earlier mistake in which I accidentally reversed the color scheme.

surfdatalegend3.mw

And here is surfdata and a legend (with a few alternatives to look & feel) for your uploaded data:

surfdataLegend_ac.mw

If you want to add up a finite number of numeric quantities like that then you are often better off using add than sum.

With beta being a Vector you cannot evaluate the indexed reference beta[k] where k is (as yet) an unassigned name.

The sum command follows Maple's usual evaluation rules, in which arguments passed to procedures (here, sum) are evaluated up front. In contrast, the add command has special evaluation rules, which allow it to defer evaluation of the first argument until such time as the index name k takes on actual values.

RBF_Interpolation_ac.mw

If you do load not either the Units:-Simple or the Units:-Standard packages then the usual operators like `+`, `*`, etc, will not combine or merge units, or convert to base units (default is the SI system).

But even in that case that you have loaded neither of those systems you can still merge units through suitable use of either the combine or simplify commands.

restart;

with(Units:-Simple):

V1 := <6*Unit(mA)>:
V2 := <3*Unit(V)>:

res2 := V1 *~ V2;

                [ 9         ]
                [--- Unit(W)]
                [500        ]

restart;

with(Units:-Standard):

V1 := <6*Unit(mA)>:
V2 := <3*Unit(V)>:

V1 *~ V2;

                [ 9         ]
                [--- Unit(W)]
                [500        ]

restart;

V1 := <6*Unit(mA)>:
V2 := <3*Unit(V)>:

res1 := V1 *~ V2;

          res1 := [18 Unit(mA) Unit(V)]

combine(res1, units);

                [ 9         ]
                [--- Unit(W)]
                [500        ]

simplify(res1);

                [ 9         ]
                [--- Unit(W)]
                [500        ]

units_comb.mw

What's wrong with the call to solve that you mentioned?

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):

tsols := [ solve(f(t)=g(t), t) ];

                tsols := [0, -1]

seq( [f(t),g(t)], t=tsols );

      [2, 2], [exp(-1) + exp(1), exp(-1) + exp(1)]