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

Because of operator precedence rules (see ?operators,precedence), you can't replace assuming with &as in such a way that the usage will be completely identical. The & operators have a very high precedence. But it'll work if you use parentheses and brackets, like this:

`&as`:= eval(`assuming`):
[is(a^2 > 0)] &as (a>0);

     true

So the first argument must be enclosed is square brackets, and the second argument in parentheses.

The RootOf form is done like this:

alias(sqrt2= RootOf(x^2-2));

And then you use sqrt2 wherever you would've used sqrt(2). There are many examples of this at ?factor, ?evala, and most of the subpages of ?evala. Whereas this RootOf form is required by many of the evala commands, the help that you quoted means that the PolynomialIdeals commands will also work with the simple sqrt(2).

In this Answer, I will assume that you can figure out how to import your matrix from the text file into Maple (for example, by using the ImportMatrix command). If you have trouble with that, just ask. And I will assume that the matrix is named A. Then do,

plots:-surfdata(
     Matrix((21,21), convert(A[..,3], list)),
     0..10, 0..10,
     labels= [X,Y,``], axes= box
);

The 21s are because you have 21 values of X and 21 values of Y. The 3 is to select the 3rd column of the matrix. The 0..10s are the ranges of X and Y.

If your Y is actually between 0 and 20 with a step of 0.5, then change the second 21 to 41, and change the second 10 to 20.

The parameters of your pdf are d, b, and r. Given numeric values for two of the parameters, we can solve (numerically only) for the third by finding the value that makes the sum for all x equal to 1. In the code below, I chose to set b and r and solve for d. Then I Sample the resulting pdf.

restart:
B:= (d,b,r,x)->
     binomial(x+r-1, x)*((d*x+1)/(1+d*x+(1/2)*b))^x*((1/2)*b/(1+d*x+(1/2)*b))^r/(d*x+1):
SumB:= proc(d,b,r) local x; evalf(Sum(B(d,b,r,x), x= 0..infinity)) end proc:
solve_for_d:= (b,r)-> fsolve(d-> SumB(d,b,r) - 1, 1.):
(b,r):= (2,2):
d:= solve_for_d(b,r);

B1:= Statistics:-Distribution(
     Type= discrete,
     ProbabilityFunction= unapply(B(d,b,r,x), x),
     Support= 0..infinity,
     DiscreteValueMap= (n-> n)
):
S:= Statistics:-Sample(B1, 50, method= [discrete, range= 0..100]):
convert(S, list);

@Josolumoh Sorry, I was reacting to what I perceived to be an unnecessary re-asking of your main request

Please, how can I do simulation study around the function?

combined with what I perceived to be a whiny "Please". I was probably overreacting. Yes, you've answered the questions promptly.

Anyway, I will show you an example of sampling a discrete custom distribution. I won't use your pdf because of the problem with the summation that Markiyan just pointed out. Instead, I will pretend that Maple doesn't already know the binomial distribution and I will create it from scratch.

restart:
B:= x-> binomial(n,x)*p^x*(1-p)^(n-x):
B1:= Statistics:-Distribution(
     Type= discrete,
     ProbabilityFunction= B,
     Support= 0..n,
     DiscreteValueMap= (n-> n)
):
(n,p):= (23,1/3):
Statistics:-Sample(B1, 10, method= [discrete, range= 0..n]);

If you get the bugs worked out of your distribution such that the pdf sums to 1, then I'll redo the example using your function.

The solve command returns two solutions for any k. Your printf statement can't handle complex numbers. Don't blame solve for that. If you want to see solve's output, use print rather than printf. Also, I don't know what you think eval is doing in your code, but in reality it isn't doing anything. Perhaps you meant evalf?