If you were looking for all non-isomorphic graphs with various properties, there is a helpful command GraphTheory:-NonIsomorphicGraphs but unfortunately for you, there doesn't seem to be any machinery which deals specifically with tree isomorphism (and clearly two non-isomorphic trees might be isomorphic as graphs).  That said, it shouldn't be too difficult to build a procedure to build them all recursively.

Something like this:

acc := proc(a, bo)
local b;
    b := bo;
    if b=3 then
       b := 2;
    end if;
end proc;

Look at the help page ?type/Matrix a lot of what you are trying to do should work, but not the restrictions on the dimensions.  You will have to check the dimensions of the inputs in the procedure body.

The error is coming from sum() which does formal summation rather than straightforward addition.  Your first step should be to replace sum() with add().

The next fix is to just do away with f_list,

eval(f,x=(0+((200-0)/100)*j))

does the same thing as your eval of f_list (which has an indexing error in it as written in the OP).

In some cases, you might be able to fake the inequalities in fsolve by using bounds on the variables:

fsolve( expr, x=0..value);

The other option is to the the feasible point finding facility in Optimization using a constant objective function:
 

Optimization:-Maximize(1, { eqns, ineqns});

There is a limitation here, however, that you cannot use strict inequalities.

There is a bug in 2D input that makes shift+enter work in unexpected ways inside lists and other enclosed enivronments.  My work around for this is the go remove the closing bracket in a list or paren in a function call (etc), then shift+enter works as expect.

That is, the first example below does not work, and second does:

[1,<shift-enter> 2]
[1, <shift-enter> 2

In this case, solve can confirm:

solve({r > 0, s > 0, r^2-r*s+s^2 > 0});

The solve command here uses SolveTools:-SemiAlgebraic to solve/simplify the set of inequalities.

Lists in Maple are immutable datastructures.  Any changes you make to a list are illusory and involve creating a new list.  Try it!
 

A := [1, 2, 3];
B := A;
A[1] := 4;
print(B);

vs.

A := Array([1, 2, 3]);
B := A;
A[1] := 4;
B;

The Data Structure sections of the Maple Programming Guide are essential reading on this topic.  Understanding what is going on here is pretty crucial to being able to write high performance in Maple.

You can gointo the View menu and turn off the autoexpansion of sections on execution.  By default, it is turned on (shown).

You might find this handy:
 

r := dismantle['string'](x+sin(x));

You can efficiently generate the set of equations using the seq command.

n := 100;
h := 1/(n+1);

# boundary conditions
bdc := {seq([u[0, j] = 0, u[n+1, j] = 0, u[j, 0] = 0, u[j, n+1] = 0][], j = 0 .. n+1)};

# interior points:
eqs := eval(
{seq(seq(
(u[i+1, j]-u[i, j])*(u[i+1, j]-2*u[i, j]+u[i-1, j])+h*(u[i,j+1]-2*u[i, j]+u[i, j-1]) = 0,
 j = 1 .. n), i = 1 .. n)},
 bdc);

However, this is a very large system of quadratic equations and probably isn't condusive to having a general purpose solver applied to it.  fsolve, for example, can solve this for n=6, but gives up after 90 seconds without finding an answer when n=7 (49 equations).  At n=100 you have 10000 equations and probably no hope with fsolve returning at all.

It looks like you might be able to solve your equations symbolically in terms of C__A

a := tau = (C__A0-C__A)/(-r__A);
b := tau = (C__B0-C__B)/(-r__B);
c := tau = (C__C0-C__C)/(-r__C);
d := tau = (C__D0-C__D)/(-r__D);
sol := solve([a, b, c, d], {C__B, C__C, C__D, tau})

I can't tell for sure since you don't provide the definitions for all the other symbols not solved for.
Even if you can't solve symbolically you might be able to wrap the call to fsolve in a procedure and plot that procedure.

@gkokovidis's soltuion does most of the job, but assume, Digits, and assign are not needed here, and the choice of i=0..3 does not need to be made appriori, if you use isolve:

eqns := {10*cos((6*(1/10))*t)-10*cos((3/10)*t+(1/4)*Pi),
    10*sin((6*(1/10))*t)-10*sin((3/10)*t+(1/4)*Pi)};
tsols := solve(eqns, allsolutions);
zsols := isolve({rhs(tsols[]) > 0, rhs(tsols[]) < 70});

sols := {seq(eval(tsols, s), s in zsols)};

sols in this case:

The command you want is probably plots:-surfdata you can then add a rotation annimation with the plot3d/viewpoint option.

If your data is in a Matrix M, something like this should work:

plots:-surfdata(M, viewpoint="circleleft");

Assignment and evaluation is a very fragile way to specialize the variables in expressions - and it is especially bad when paired with the wierd evaluation rules for Arrays/Matrix.  I prefer to use eval, which works well here:

M__11 := Matrix(eval(Mm__11,
    {I__0 = 0.102033222998e-2,
     I__1 = .80922342232909,
     I__2 = 1,
     a = .2223212330091,
     b = .293939202032292
    }), datatype = float[8]);