Fern := proc(n::posint
             , x::Array(datatype=float[8])
             , y::Array(datatype=float[8])
             , r::Array(datatype=float[8])
            )
local i, pick;
option autocompile;
    for i to n-1 do
        pick := r[i];
        if pick < 0.01 then
            x[i+1] := 0;
            y[i+1] := 0.16*y[i];
        elif pick < 0.08 then
            x[i+1] := 0.2*x[i] - 0.26*y[i];
            y[i+1] := 0.23*x[i] + 0.22*y[i] + 1.6;
        elif pick < 0.15 then
            x[i+1] := -0.15*x[i] + 0.28*y[i];
            y[i+1] := 0.26*x[i] + 0.24*y[i] + 0.44;
        else
            x[i+1] := 0.85*x[i] + 0.04*y[i];
            y[i+1] := -0.04*x[i] + 0.85*y[i] + 1.6;
        end if;
    end do;
    NULL;
end proc:

n := 10^5:
r := rtable(1..n, frandom(0..1,1), 'datatype'=float[8]):
x := Array(1..n, 'datatype'=float[8]):
y := Array(1..n, 'datatype'=float[8]):

(Fern(n,x,y,r));
plot(x,y, 'style=point', 'color=green');

Using assign prevents you from reusing the assigned variables in the set of equations.  A better approach is to symbolically solve the set of equations, then create a procedure from that solution, and plug each random combination into it.  For example,

restart;
with(combinat):

vars1 := {a, b, c, d, e, f, g, i, m}:
vars2 := [h, j, k, l, n, o, p, z]:
sol := solve({NULL
              , a+b+c+d-z = 0
              , e+f+g+h = z
              , i+j+k+l = z
              , m+n+o+p = z
              , a+e+i+m = z
              , b+f+j+n = z
              , c+g+k+o = z
              , a+f+k+p = z
              , d+g+j+m = z
             }
             , vars1):

f := unapply(map(rhs,sol) union {vars2[]}, vars2);

for cnt do
    A := randcomb(100,8);
    if nops(f(A[])) = 16 then break;
    end if;
end do:

cnt; f(A[]);

It isn't clear what you really want.  What is the expression you want to plot?  The triple integral you show fully evaluates to a number. 

I didn't see a mention of the RationalMap procedure in that paper.  Did you mean ?RationalMapImage or ?RationalMapPreImage, which exist in the ?RegularChains[ConstructibleSetTools] subpackage.

You should write this as a procedure with a single argument, the Matrix to check.  The procedure should return either true or false. Here's a straightforward implementation, with enough comments to explain what is happening:

IsDiagonallyDominant := proc(M :: Matrix)
local strict,sumrest,diag,i,j,m,n;
    strict := false;   # flag to keep track of whether at least one row is strict
    (m,n) := op(1,M);  # extract number of rows (m) and columns (n) of Matrix M
    for i to m do
        diag := abs(M[i,i]);  # Get the magnitude of the diagonal term
        # Sum the magnitudes of the other terms in the row.
        sumrest := add(abs(M[i,j]), j=1..i-1) + add(abs(M[i,j]), j=i+1..n);
        if diag < sumrest then
            return false;  # short cut; return false if row fails test
        end if;
        # Assign strict, if needed
        if not strict and sumrest < diag then
            strict := true
        end if;
    end do;
    # At this point, none of the rows fail.  
    # To be diagonally-dominant, at least one row 
    # must meet the strict equality,
    # so return the value of strict.
    return strict;
end proc:

Here's a somewhat interesting extension.  Rather than return just true or false, and fail if an element is non-numeric (or, say, Pi, which is a deficiency in the above routine), we can extend this to return a boolean expression:

IsDiagonallyDominant := proc(M :: Matrix)
local pred,strict,sumrest,diag,i,j,m,n;
    pred := true;
    strict := false;
    (m,n) := op(1,M);
    for i to m do
        diag := abs(M[i,i]);
        sumrest := add(abs(M[i,j]), j=1..i-1) + add(abs(M[i,j]), j=i+1..n);
        if diag :: realcons and sumrest :: realcons then
            if evalf(diag < sumrest) then
                return false;
            end if;
            if evalf(sumrest < diag) then
                strict := true;
            end if;
        else
            pred := pred and sumrest <= diag;
            if strict <> true then
                strict := strict or sumrest < diag;
            end if;
        end if;
    end do;
    return pred and strict;
end proc:

M := <<2|1>,<c|d>>:
IsDiagonallyDominant(M);
                                       | c | <= | d |

That is not a well-formed Maple expression.  The = and > operators have the same precedence and  are non-associative (see ?operator,precedence), so you need to add parentheses to group the components. Regardless, I don't see how false could be returned. 

How did you define the functions f[i]?  Are they really Maple procedures, for example,

f[1] := x -> x^2:

Note that E[ij] should be E[i,j]. 

myproc := proc(x0)
local i,x;
  x := x0;
  for i from 0 to 10 do
      x := x-(some stuff);
      print(x);
  end do;
end proc:

How 'bout plotting a mold for it:

cylx := y^2+z^2-1:
cyly := z^2+x^2-1:
cylz := x^2+y^2-1:

expr := cylx*cyly*cylz;
plots:-implicitplot3d(expr, x=-1..1, y=-1..1, z=-1..1, numpoints=10\000);

A minimal example that demonstrates the problem would be helpful.

The problem is that you are attempting to assign to a parameter. That generally is not allowed in Maple.  What do you hope to accomplish with this procedure?  If you merely want to assign to the global variable x, then you could do

testproc(x0,y0)
global x;
    x := y0;
end proc:

While you could use ?assigned to determine whether a name is assigned, you might be better off using ?PolynomialTools[CoefficientList] and directly combining multiple calls.  Say

polys := { a2*x^2 + a1, b2*x^2+b1, c1*x }:
eqs := {seq(PolynomialTools:-CoefficientList(p,x)[], p in polys)};

You can use ?expand to handle this.  However, calling expand directly does not do what you want:

y := cos(2*x+(Pi/6)):
expand(y);
                                      1/2
                       1/2       2   3
                      3    cos(x)  - ---- - sin(x) cos(x)
                                      2

The problem is that expand was too aggressive.  To dampen its enthusiasm, you can pass it optional arguments that tell it what not to expand:

expand(y, 2*x);
                                     1/2
                       1/2 cos(2 x) 3    - 1/2 sin(2 x)

The reason sum does not work here is that the expression being summed, specifically the derivative, is evaluated before the index has been assigned a value.  You could use sum, though add is the better choice, if you delay the evaluation of the diff function.  That can be achieved with forward quotes:

sum('diff'(f, B[j])^2, j = 1 .. 2);
                                     2    2
                                    b  + a

If the signal is in degrees, then divide by 0.3, use the Real-To-Integer block, then scale the output by 0.3.  Generally an angle will be in radians, so you will want to adjust the scale factors accordingly.