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?

Here's an easy way that I think corresponds to what you were asking for. We sort the points returned by LagrangeMultipliers based on the function value, lowest to highest. Then we pair each point with its function value. This requires Maple 18 or later. If you have an earlier Maple, let me know---it can still be made to work.

F:= (x,y,z)-> x*y*z:
Sols:= Student:-MultivariateCalculus:-LagrangeMultipliers(
     F(x,y,z), [x^2+4*y^2+4*z^2-4], [x,y,z]
):
Sols:= sort([Sols], key= evalf@F@op):
map(s-> [F(s[]), s], Sols);

Any data structure can be stored in a Vector. But I think that in your situation storing them in a list would work better. I'm not sure about your situation, and I'd like to see more details.

Here's a procedure that does the job with simply three sorts and a reindexing. It doesn't use any conversion to disjoint cycle form; it works strictly with the permlist form.

Perm1to2_D:= (L1::list, L2::list)->
     sort(L1, 'output'= 'permutation')
          [sort(sort(L2, 'output'= 'permutation'), 'output'= 'permutation')]
:

Here's how it works. Given a permutation in permlist form, its inverse can be obtained by sorting it with output= permutation. Given two permutations P1 and P2 in permlist form, the product or composition of P2 applied to P1 can be obtained by simply indexing: P1[P2]. 

Here's a time test of the four procedures so far presented for this problem.

Download permL1toL2.mw

The first argument to animate shouldn't be a plotting command; rather, it should be just the procedure name of a plotting command. The second argument should be a list of arguments for the procedure in the first argument parametrized by the animation parameter. The third argument should be the animation parameter and its range. So, your command should be

plots:-animate(
     plots:-dualaxisplot,
     [sin(x+A), cos(x+A), x=0..5, style = line, gridlines = false],
     A = 0 .. 5
);