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.

The conversion of an expression sequence to a list and vice versa is a triviality both in terms of coding and efficiency. Kitonum has shown the coding syntax. Regarding the efficiency, both the time and memory needed for the conversions are very small constants; they do not depend on the length of the data.

It is often better to stick with lists because of the indexing problems that can occur when an expression sequence has 0 or 1 elements.

This should not have been a separate Question, but since Tom already put an Answer here, I won't delete it.

To me, the input diff(x-1 = 0, x) is not GI (garbage input) nor is 1=0 GO (garbage output). As I was trying to explain in the other thread, everything makes sense if you take into account the quantifications that are implicitly present.

Maple allows for the simultaneous differentiation of both sides of an equation, and this is often useful. However, you could argue that this shouldn't be allowed (because of this very problem). And thus you could argue that diff(a = b, x) should be considered GI and an error for any a and b. If this restriction were to be imposed, I wouldn't object. It's trivial to use overload to impose that restriction:

restart:
diff_orig:= eval(diff):
unprotect(diff):
diff:= overload([
    proc(e::`=`)
    option overload; 
        error "Differentiation of equations not allowed."
    end proc, 
    diff_orig
]):
protect(diff, diff_orig):

The output 1=0 could be viewed as an algebraic equation or as a boolean relation. If it's an equation, then it's an equation whose solution set is empty. That's not garbage. If it's a relation, then it evaluates to false; false is not garbage.

The input and output considered together constitute a proof-by-contradiction that there does not exist any differentiable function f such that f(x) = x - 1 and f(x) = 0 (for all x in whatever domain of differentiation you want to consider). That's knowledge (albeit trivial), and knowledge isn't garbage.

You could do it in the way that you suggest if you want, but all the standard differential operators are in the VectorCalculus package: Gradient, Jacobian, Hessian, and several others.

As a workaround, you can do this:

with(Physics):
Physics:-Setup(assumingusesAssume= false):
csgn(a - b[2]) assuming a > b[2];

The simplify and all the other assumptions are superfluous. You can put them if you want, but they make no difference.

The number of combinations that you're trying to generate is so large that they could not be listed by any computer (or any other method), ever. And it's not possible that there ever will exist a computer (or any other method) of doing this. We can count the number of such combinations like this:

Number(Iterator:-MultiCombination([2^16 $ 2^8], 2^16));
7820701318841130906613553459164687216584467046464667050034224507
  41484858519585179827709391651098624501404918769795637653274401
  57029711677787977640862663017116097938315215722518761350936266
  38516727105642689147056759932239793313894413127816886156554817
  76804366751588302164524267450481108817377756113734373353127476
  60436555681987839580644837879373377115558668241858194605420318
  57090864326621621084868745055201958236210246318706997513184396
  62089540815829177821487234240936427010441020860670200642172134
  59286677922608980438360909995794377076419365877451789293653200
  24753539064997228881393914518695160243181091267722013448011707
  00847339475028997519306136404337673374805865549565962012504690
  7222757907787394934267401909783755764225

evalf(%);
                                     723
                       7.820701319 10   

So, if every cubic meter in the entire Universe contained a computer, and they each generated one combination each nanosecond over the entire age of the Universe, then the progress they would've made on the job would still be far less than the progress made on the job of counting every atom in the Universe by counting a single atom.

So, you need to say what you want to accomplish by generating these and figure out a way to do it without generating them.

The way that you're using indexing (numbers appearing in square brackets appended to expressions) makes no sense to me, and the error message implies that it makes no sense to pdsolve either. I have no idea whether pdsolve will ultimately be able to solve your problem, but at this point it doesn't even recognize the problem as a PDE.

Suggestions:

1. Don't use 2D input.
2. Make your primary goal to get pdsolve to recognize your system as a PDE without error. Having the system prettyprinted on screen in a certain way should only be a distant secondary goal.
3. Don't use PDEtools:-declare because its purpose is primarily to address cosmetic rather than syntactic issues.
4. If you want 2 us (or any other variable), just call them u0 and u1. 
5. Stop saying that you want to solve it by the perturbation method (or any other method). If you can't even input the problem, then the method chosen to solve it is at the moment irrelevant.

The cosmetic issues can be addressed afterwards. They are of course irrelevant if it can't even recognize what you've entered as a PDE.