Example:

P1:= plot3d(x*y, x= -1..1, y= -1..1):
P2:= plots:-pointplot3d([[0,0,0]], symbol= solidsphere, symbolsize= 32, color= red):
plots:-display([P1,P2]);

Why is your data organized in this strange way? Something tells me that you're making it more complicated than it needs to be. If each letter is to be paired with a number, then store them as pairs:

map(L-> {seq(L[2*k-1]=L[2*k], k= 1..nops(L)/2)}, data);

You see that when you use this more natural format, the sets are automatically sorted by letter. Unlike the other two Answers, mine will work on sublists of any length.

Your manual solution is faulty because the Laplace parameter s still appears in the equation. Since no gamma appears in your solution, I will assume that the coefficient gamma in the second original equation was supposed to be lambda.

These can be solved easily with Maple, like this:

restart:
Ps:= [P[0],P[1]]:
eqs:= Ps(s) =~ [1/(s+lambda)+mu*P[1](s)/(s+lambda), lambda*P[0](s)/(s+mu)]:
Ls:= solve(eqs, Ps(s))[];

Ps(t) =~ inttrans[invlaplace]~(rhs~(Ls), s, t);

If the step size is (b-a)/n, then including a and b, there are n+1, not n, evenly spaced points in the interval.

They can be generated like this (I take the name linspace from the equivalent Matlab function):

linspace:= (a,b,n)-> [seq(a+(b-a)*k/(n-1), k= 0..n-1)]:
linspace(0,1,50);

The procedure f can be applied to them like this:

f~(linspace(0,1,50));

While time is useful for testing a block of code (which should be in a single execution group), it is better to use CodeTools:-Usage to test a single procedure call. It returns more-detailed information. Example:

P:= `*`('randpoly(x, degree= 2^9, dense)' $ 2^7):
Prod:= CodeTools:-Usage(degree(expand(P)));
memory used=0.58GiB, alloc change=257.00MiB, cpu time=16.59s, real time=16.62s, gc time=0ns

     Prod:= 65535

By default, Usage returns directly the result of the computation and returns its usage report as a side effect. This default can be changed with Usage's output option.

This is a common source of bugs in Maple code: The code is expecting an expression sequence. If the sequence has only one member, then indexing the sequence becomes indexing of the member, which is not what was intended. The solution is to turn the expression sequence into a list by enclosing it in square brackets:

Answers:= [solve(equation that needs solving)];

This will work even if solve returns nothing. Now if Answers has only one solution, Answers[1] will return that solution in its entirety.

Terminology: Neither a list nor a sequence should be called a matrix, even when it can be double-indexed. And an expression sequence shouldn't be called a list, as in the title to this Question.

The number of elements of a list is returned by nops. This'll tell you what are the safe indices to use.

You should not use indexing to extract the a and b values from an individual solution. Use eval instead, as in

eval(a, Answers[1]), eval(b, Answers[1]);

Maple has tools to make the elimination of duplicates trivial and the elimination of near duplicates very easy.

Your code in both this Question and your previous Question seems to suggest that you have the mistaken belief that a data structure's indices need to be integers. This isn't true for a Maple table, which is what your x is. Its indices can be anything at all; they don't even need to be numbers. In the code below, I use the roots themselves as the indices, which automatically removes exact duplicates. Your attempt to use the integer indices i and j makes your code more complex.

Download HankelRoots.mw

You have at least two errors with your final integral. One is a logic error, and the other is a syntax error. The logic error is that you use both x and y as variables and that you integrate with respect to x. I don't understand why you introduced x at all, and certainly the integration should be performed with respect to y.

The syntax error is that ay is being treated as a single variable. If you intend multiplication (as I'm sure you do), then you need to enter either a space between them (a y) or an explicit multiplication operator (a*y).

The issue with the decimal points is something called floating-point contagion: Once there's one decimal point in an expression, it tends to spread to other numbers in the expression automatically as operations are performed. In your case, the one decimal point comes from the 62.4. There are two ways you can prevent floating-point contagion. The first (and my preference) is to represent floating-point constants as symbols; the second is to represent them as exact rationals such as 624/10.

Use the seq command to create the lists sys and var, like this:

N:= 10:
sys:= [seq(galerkin_funcs[i], i= 1..N+4)]:
var:= [seq(w[i], i= 1..N)]:
(Kmat, Fmat):= LinearAlgebra:-GenerateMatrix(sys, var)

I did this without a ruler, just eyeballing it and restricting my coordinates to integers. You should be able to see how to make any numeric adjustments.

rectangle:= (P::[numeric,numeric], h::numeric, w::numeric)->
     plots:-polygonplot([P, P+~[w,0], P+~[w,h], P+~[0,h], P], args[4..])
:
     
plots:-display(
     [plot([[0,0], [10,0], [10,10], [0,10], [0,0]], thickness= 5, color= black),
      plots:-polygonplot([[0,10], [5,15], [10,10], [0,10]], color= red),
      rectangle([1,1], 3, 4, thickness= 2, color= white),
      rectangle([7,0], 4, 2, thickness= 2, color= blue),
      rectangle([1,7], 2, 2, thickness= 2, color= pink),
      rectangle([7,7], 2, 2, thickness= 2, color= pink),
      plots:-polygonplot([[8,12], [8,13], [9,13], [9,11], [8,12]], thickness= 2, color= grey)
     ], scaling= constrained, axes= none    

);

The following appliable module will extract the minimal and maximal points (as determined by the final coordinate) from any two- or three-dimensional plot. If an extremum is achieved at multiple locations, only one location is reported.

PlotExtrema:= module()
local
     Min, Max, MinPt, MaxPt,
     ModuleApply:= proc(P::specfunc(anything, {PLOT,PLOT3D}));
          Min:= infinity; Max:= -infinity; MinPt:= []; MaxPt:= [];
          plottools:-transform(`if`(op(0,P)=PLOT, F2D, F3D))(P);
          MinPt, MaxPt
     end proc,
     F2D:= proc(x,y) Scan(args); [x,y] end proc,
     F3D:= proc(x,y,z) Scan(args) end proc,
     Scan:= proc()
     local Pt:= [args];
          if Pt::list(numeric) then
               if Pt[-1] < Min then
                    Min:= Pt[-1]; MinPt:= Pt
               end if;
               if Pt[-1] > Max then
                    Max:= Pt[-1]; MaxPt:= Pt
               end if
          end if;
          Pt
     end proc
;
end module:

Examples of use:

P:= plot3d(x*y, x= -1..1, y= -1..1):
PlotExtrema(P);

     [-1., 1.00000000000000, -.999999999999999], [-1., -1., 1.]

P:= plot(x^2+1, x= -2..1):
PlotExtrema(P);

     [0.514891618090418e-2, 1.00002651133784], [-2., 5.]

Of course, my algorithm only uses points that are actually plotted!

Your resulting equation is wrong: It should be x+y = z. You can do the division in one command as

normal(eq1 / eq2);

The division operation `/` "zips" over equations. The normal simplifies the quotients (so that (x^2-y^2)/(x-y) becomes x+y).

To see the code of `factor/factor`, do

showstat(`factor/factor`);

Note the extra quotes, which are required for any name that contains a special character, such as /. That's probably most internal procedure names.

How type(..., rational) works is trivial, there's no sublety about it: An object is of type rational if and only if it is explicity and manifestly an integer or a ratio of integers. So, if I construct some radical expression expr which can be simplified to an integer, but doesn't simplify automatically, then type(expr, rational) would return false. So there is definitely no reason why you "should believe that type(x, rational) would work when the description of x is quite complicated." Likewise, I wouldn't be surprised if you could create a cubic polynomial whose coefficients can be simplified to integers but which aren't manifestly integers.

I think that you have a completely wrong idea about what type is supposed to do. The expression type(x, T) doesn't "test whether" x "evaluates to"* an expression of type T; it checks whether x is an expression of type T. In other words, type is a syntax checker, not a mathematical-set-membership checker. If you're looking for the latter, the command is is. (Its prowess is often debated here, and it's certainly not always correct, but at least in its intention is is the mathematical-set-membership checker.)

*Before someone jumps in here to say that I'm using the word evaluates incorrectly, let me point out that I am quoting the OP's usage, which was intended in a mathematical sense of the word, not the way that the word is used by Maple.

This was covered somewhat in my Answers to your previous Question about permutations. But, for a quick review:

L:= [1,2,3,4]:
P:= [[1,3], [2,4]]:
L[convert(P, permlist, 4)];

The third argument to convert, 4 in this case, is the degree of the permutation. The command won't simply assume that the highest number that appears is the degree because, for example, 5 could've been a fixed point of the permutation.

The procedure that generates the random number must return unevaluated when passed a symbolic argument and must be declared with dsolve's option known. The error control must be limited lest the numeric solver will try to take infintesimal steps because of the chaotic behavior. One way of doing that is with dsolve's option minstep. Putting it all together, this works (I made up some arbitrary initial conditions for your equation):

eq:=diff(a(t),t,t) + a(t) = 3.*(1+0.1*R(t)/1000000000000.)*cos(5.*t):
R:= t-> `if`(t::complexcons, rand(), 'procname'(args)):
Sol:= dsolve({eq, a(0)=0, D(a)(0)=1}, numeric, known= R, minstep= .001);