@Carl Love I have no problem to use for. I was wondering if Maple had something similar to NestList (Mathematica).   

@Carl Love Many thanks.  Nice piece of code.  And if I am only interested in having the map symbolically iterated (I am avoiding using floats).  

@vv Many thanks but I am already using (f@@p) to find the fixed points.   What I need is to have the system iterated as in the initial post. More than that, I need to have it iterated backwards (two solutions each time). 

@Kitonum Many thanks again. If rsolve was not available, would it possible to implement the iterations without a for, that is, using commands such as seq?  

@Kitonum Many thanks.  It helped a lot.   And when it is necessary to iterate the map?  

@acer Many thanks.   I wasn't aware of the command Alias.  

@Carl Love Many thanks for the procedure.  As for your previous question, the example is indeed to show that round-off errors can lead to anomalies but can be also used for calculating some measures (largest positive lyapunov exponent for instance).  

@Carl Love Many thanks, Carl.   I got an error message when issuing the commands:

Hx:=`hf<->hex`(0.25);
                    Hx := "3FD0000000000000"
Hf:=`hf<->hex`(Hx);
Error, (in sscanf) expecting exactly 16 hexadecimal digits

Please note that the hexadecimal representation is different from what appears in the original message due to the float precision.  Would it be possible to simulate 16-bit and 32-bit floats?

As for the addition procedure, yes, please (I always learn a lot from the way you code the solutions).  
 

@John Fredsted OK, I go it.   Many thanks.   Following your ideas, I can build myInverse to block invert the matrices.  

@John Fredsted Many thanks.   I have just copied and pasted the example you've just send and myTranspose returned the answer for the variable expr (defined earlier) and not for M.  What am I missing?  

@John Fredsted Many thanks.   Do I always have to define the matrix to be defined as expr?   

@Carl Love Many thanks. a)  I had that and the solution figured out before you sent the reply, but it is good to know that I got it right.  b) I have completely forgotten that.  Many thanks.  

@Carl Love Dear Carl,  thank you very much (again).

If you have the time could you help me with the following issues?

a) OK.   I did as suggested and some of the results I was expecting to see showed up.  Example: eval(CDF(X2,v),[v=0.64,m=3]).  However when I tried eval(Quantile(X2,v),[v=0.95,m=3]), Maple 2017 returns FAIL.  That is a bit discouraging since I need to increase m up to 50 and calculate the quantiles. I hope I am missing something here.

b) Could you explain what does $ 2 do?   

Again many thanks.

@Carl Love Many thanks.   As you know I work with sets of polynomials.  For each of them I need to make a call to Groebner independently. If I got it right, I can send each call to a core (24 core each time) and them save the results, but unfortunately I don't know how. 

@Carl Love models=subsets. I need all the subsets of [f(x,y,z),g(x,y,z),h(x,y,z)] with f(x,y,z) containing that specific monomial. Please note that MyTest will remove some of the subsets. I don't know the package but it sounds worth trying. The functions f,g and h are full quadratic polynomials in x,y and z (V as used in the procedure). Many thanks.