@Carl Love Many thanks again.  My bad! My copy and paste wrecked your beautiful code.  I do apologize for that.  I have installed Maple to a new computer and forget to change the settings (you warned me about it years ago).  I will try not to do that again in future posts.

Thanks for the piece of code.  The only part I am not sure if when you calculate aaa.  Why 49? 

When using hf the number of digits is set to 15 and I need to have it changed accordingly (more or fewer periods).  

Why did I choose the fifth fixed point?  The first 7 values can be expressed as radicals when using allvalues.  The last five fixed points out of 7 are the values of the period-5 cycle.  I chose one of them for no particular reason. (I couldn't get Mathematica to show those values).

My biggest concern is that I used evalf to convert the initial condition, that is fp[5], to a float number (but the result could be complex with 0*I for some reason).  My guess was that Maple would use float all the way to the end and, since there are only multiplications of numbers smaller than one, there is no reason to think that Maple would come up with infinity or a complex number.   (My plan is to have only one NestList to do both symbolical and numerical iterations).  

@nm Many thanks.   I have tried unapply but I was definitely missing something. Actually I was adding stuff that shouldn't be there.  

@Carl Love Many thanks again.   

@Carl Love Thank you ever so much.  I had no problem to run them all.  Many thanks again.  

I have encapsulated the second solution as

logistic := y-> 4*y*(1-y);
logistic1:=proc(f,x,n::nonnegint)
local res,p,x0;
x0=x;
res:=seq((x0:=f(x0)),p=1..n);
return [res];
end proc:

When I issued logistic1(logistic,1/4,2), it worked.  Many thanks.

@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?