Edit: I just read your followup to Rouben that you posted while I was writing this, so I understand now that you want to do more than just integrate this function. Getting your worksheet now.

Original answer: The function is piecewise polynomial with exact rational coefficients and breakpoints. It's trivial to integrate it symbolically (even if it weren't smooth or continuous) and get the exact value.

int(Z(x), x= 5/3..5);

             361954048721879 / 361954048856400

Approximating the function by some other function is a bad idea. The computational effort required to get an approximating function is far greater than that required for integration.

You have defined

G:= array(1..N, 1..N);

You should change that to

G:= Array(1..N, 1..N);

Then G will print. The command array is deprecated; you should never use it.

Here is a much more general way by which we can work with differential operators whose behavior mimics that of D itself, and when applied to arguments, the expressions are converted to diff form. I'n going to call my operator DiffOp, but feel free to change that name.

`evalapply/DiffOp`:= proc(DO, Y)
    local y,x;
    if procname::indexed then
        y:= op(procname); x:= Y[];
        evalindets(
            evalindets(
                op(DO),
                {
                    'D'(identical(y)), 
                    And(
                        typefunc(
                            identical(y), '`@@`'(identical(D), anything)
                        ),
                        function(identical(y))
                    )
                },
                d-> convert(d(x), diff) 
            ),
            identical(y),
            f-> f(x)
        )            
    else
        y:= Y[];
        DiffOp[y](
            subs(
                {D@@0, `1`}=~ y,
                evalindets(
                    evalindets(
                        subs(D= D(y), op(DO)),
                        '`@@`'('D'(identical(y)), anything),
                        d-> (D@@op(2,d))(y)
                    ),
                    {procedure, And(appliable, `module`)},
                    p-> p(y)
                )
            )
        )
    fi
end proc
:

#Example: 
# Use D@@n instead of D^n. The @@ is Maple's operator for functional iteration. 
# D^n will be interpretted as the nth algebraic power of the first derivative.
# Use `1` or D@@0 (the choice is yours) for the identity operator (which will become
# y(x) when the differential operator is applied.

did:= x*D@@3 + 144*sin*D@@12 + 7*D + `1`;
This output makes more sense when prettyprinted in a worksheet:
        did := x @@(D, 3) + 144 sin @@(D, 12) + 7 D + 1

DiffOp(did)(y)(x);

  /  3      \                 /  12      \                      
  | d       |                 | d        |     / d      \       
x |---- y(x)| + 144 sin(y(x)) |----- y(x)| + 7 |--- y(x)| + y(x)
  |   3     |                 |   12     |     \ dx     /       
  \ dx      /                 \ dx       /                      

                 

I will explain to you why this doesn't work with eval or subs, and I'll even teach you an alternative that does work, but I must emphasize that your whole idea along these lines is very, very, VERY cruddy and prone to error; so, you must never do this or ask a Question about it again. Your subsequent Question, already Answered by Kitonum, shows that you're well on your way to having an alternative to this mess.

Why it doesn't work: The issue is order of evaluation. An expression being modified by subs or eval is evaluated before the substitutions are made. (This is simply the default order of evaluation for the vast majority of Maple commands.) So, at the time of evaluation, i has no value. That doesn't make sense to op. There are many Maple commands that are prepared to handle situations where a parameter that is usually numeric is passed symbolic (i..e., has no value). And there are many other commands, such as op, that are not prepared to handle that.

An alternative that works: Use a procedure to change the order of evaluation:

MyC:= i-> coeff(op(i,did), D, op(2,op(2,op(i,did))));
MyC(3);

Even better:

MyC:= (i::nonnegint)-> coeff(op(i,did), D, op([i,2,2], did));
MyC(3);

Note that nested op calls can be compressed by using a list as the first argument.

Use an embedded if statement rather than `if`, like this:

str:= "A";
x:= 10;
str:= sprintf(
    "%s it was%s 10", str, if x++=10 then "" else x:= 9; " not" fi
); 

The benefit of an embedded if statement is that any number of statements can appear in its then and else clauses. It's only the value of the final of these statements that becomes the return value of the embedded statement.

Your computation can be done by very simple arithmetic in compiled code. For the n=8 case, that avoids (2^8)^3 ~ 16 million calls to Bits:-Split, (2^8)^2 + (2^8)^3 calls to ListTools:-Occurences, and the creation and garbage collection of the same number of lists. My code below does the n=8 case in 17 seconds total on one processor. Further gains can be obtained by

• Writing the code directly in C so that you can use whole-word bit operators.
• This code can be parallelized with Grid or Threads (Threads is usually faster),
• Perhaps some speedup can be achieved by using 2-byte integers directly in C.

I made a great number of changes with the arithmetic. Please verify the accuracy of this by comparing against known-correct code using the same random vector M.

restart
:
BuildA_M:= proc(
    A::Array(datatype= integer[4]), M::Array(datatype= integer[4]),
    r::integer[4]
)::integer[4];
local 
    i0::integer[4], i::integer[4], j0::integer[4], j::integer[4],
    m::integer[4], t0::integer[4], t::integer[4],
    p::integer[4], S::integer[4]
;
    for i0 from 0 to r-1 do
        for j0 from 0 to r-1 do
            S:= 0;
            for t0 from 0 to r-1 do
                p:= 1; i:= i0; j:= j0; t:= t0; m:= M[t0+1];
                while i*t <> -j*m do            
                    if (i*t - j*m) mod 2 <> 0 then p:= 1-p fi;
                    i:= i/2; t:= t/2; m:= m/2; j:= j/2
                od;
                S:= S+p
            od;
            A[i0+1, j0+1]:= S
        od
    od;
    0
end proc
:
BuildA:= Compiler:-Compile(BuildA_M):

LArip:= proc(n::posint)
local 
    r:= 2^n, 
    M:= Array(1..r, combinat:-randperm(r) -~ 1, datatype= integer[4]),
    A:= Array((1..r)$2, 0, datatype= integer[4])
;
    BuildA(A, M, r);
    (A,M)
end proc 
:
CodeTools:-Usage(LArip(8));
memory used=263.30KiB, alloc change=0 bytes, 
cpu time=16.55s, real time=16.68s, gc time=0ns

​​​​​​

The Float(undefined) could be the result of dividing by very small numbers. If you increase the precision with Digits:= 25 the problem goes away. I just guessed that 25 might work, and it did. Some smaller value may work also.

Change the line where tab_err is created to this:

tab_err:= Vector(N-1); 

Note that your code had = not :=, so it wasn't creating anything.

Change all references of tab_err[nrow, 1] to tab_err[nrow]. (There are only two such references.)

Change the plotting command to

print(plot(<<seq(1..N-1)> | tab_err>, axis[2]= [mode= log])); 

In Maple, all plotting commands do nothing more than create a data structure that represents a plot; they do not display the plot (this even includes the command named display). Displaying the plot requires a print command. The reason that plots get displayed by top-level commands is (essentially) that you automatically get a "free" print at the end of every execution group.

Almost always, a plot of the absolute values of errors vs number of iterations should be done logarithmically.

Replace the entire for loop that prints tab_err with the single command

printf("Printing tab_err\n%6.4e\nDone printing tab_err\n\n", tab_err)

Then you will see for certain that tab_err is being printed.

Before executing the procedure, give the command

Digits:= 15; #or a larger integer

The default value of 10 digits is not enough to appreciate what's happening in your procedure.

I think that what you've done is a brilliant hack that gets around Histogram's unnecessary refusal to handle the legend option.

Much of your code is redundant. The entire thing can be replaced by

plots:-display(
    evalindets~(
        Statistics:-Histogram~(
             Statistics:-Sample~(Normal~([0,1], 1), 1000),
             color=~ [gold, blue]
        ),
        specfunc(COLOUR),
        1@@0, #multi-argument identity function
        LEGEND~(["N(0,1)", "N(1,1)"])
    ),
    transparency= 0.3
);

Unlike my previous Answer to you, this one does produce exactly the same results as yours.

As shown by Kitonum, a simple use of solve proves that g must be a square except *possibly* if b*d*f = 0. So let's consider that exception. 

eq:= sqrt(a+b*sqrt(c+d*sqrt(e +f*sqrt(g)))) = h;
eq||(b,d,f):= eval~(eq, [b,d,f]=~0)[];
                               (1/2)    
                       eqb := a      = h

                                      (1/2)    
                        /       (1/2)\         
                 eqd := \a + b c     /      = h

                                            (1/2)    
                 /                    (1/2)\         
                 |      /       (1/2)\     |         
          eqf := \a + b \c + d e     /     /      = h

None of those equations contains g. If you find solutions for those (and there are many), you could use any values of g and the other missing variables. It'd be a total waste of computational effort to loop through the variables that do not appear. It's also a waste to loop through h: If the left side is an integer, then it is h, and it's the only possible h.

The best way to do this is to build the integral into the ODE system. In the following code t(s) represents your integral, so diff(t(s), s) is the integrand and t(0) = tinit is the initial condition.

#I just made up these parameter values:
g:= 9.81:  h:= 0.1:  rinit:= 0.1: tinit:= 5:

dsysA:= {
    diff(r(s), s) = sqrt((g*r(s) + h)^2 - 1 + 1/r(s)),
    diff(t(s), s) = (g*r(s) + h)/(1 - 1/r(s)),
    r(0) = rinit, t(0) = tinit
}:
dsnA:= dsolve(dsysA, numeric, abserr= 1e-6):

dsnA(0.1);
              
  [s = 0.1, t(s) = 4.89144107647636, r(s) = 0.429604740501696]

You need to be able to distiguish the free variable eps (which could even be an expression) from the variable r, which must be a variable because it'll temporarily be bound by int before being returned free (free and bound are opposites). There is no "clean" way to do this without passing r to f. By "clean", I mean without subverting the normal evaluation rules and without globally designating r as the temporarily bound variable. So, I'd do this:

f:= (g, r::name)-> int(g(r), r):

Then you could achieve your two desired calls as 

f(h, r);
f(h@(x-> x+eps), r);

This works whether g or h are defined procedures, just names, or some concoction in between those two. The shift function x-> x+eps could be predefined:

sh:= eps-> x-> x+eps:
f(h@sh(eps), r);

I recommend that you never use the syntax f(...):= ... to define a procedure (sometimes called a "function", but that word is ambiguous when talking about Maple). That syntax works sometimes, but not always.
 

If you want some numbers to be ignored, make them names rather than strings. You can do this by replacing "2222" with `2222` and likewise for the other string.

Thus, I get

By the way, your entire for loop can be replaced by map[inplace](evalf@log[10], B).

I agree that that is a bug. I've seen it many times before. It seems to be related to large fractional exponents. They cause solve to create RootOf expressions of polynomial degree the numerator of the fraction. So, in this case, that's a polynomial of degree 131537.

A workaround is to convert the exponent to a symbol:

Eq:= 1 - 1/(1+203808*exp(-342569/506*t)*(1/131537))^(131537/203808) =
44983/56599;
Eq2:= evalindets(
   Eq, anything^fraction, e-> op(1,e)^convert(op(2,e), name)
);
sol:= {solve}(Eq2, t);
evalindets(sol, name, parse);
evalf(%);

Acer's Answer overrideoption is the best way to do this. I just want to show you a primitive alternative to your primitive method.

#Your way:
c  := op(select(has, op(1, b1), COLOUR));
b2 := eval(b1, c=COLOUR(RGB, 1, 0, 0)):

#Replace with:
b2:= evalindets(b1, specfunc(COLOUR), c-> COLOUR(RGB, 1, 0, 0));

I mention this because, IMO, evalindets is the most important command in all of computer algebra and everyone should learn it.

My code does not necessarily do exactly the same thing as yours. Mine replaces all colors.