You haven't supplied any inital conditions. So, is this what you want?

DEtools:-dfieldplot(
   [diff(x1(t),t) = 3*x1(t) + 2*x2(t),
    diff(x2(t),t) = 5*x1(t) + 1*x2(t)
   ],
   [x1,x2](t),
   t= 0..2, x1= -Pi..Pi, x2= -Pi..Pi
);

Increasing the implicitplot option gridrefine from its default 1 is more sophisticated than using a large value of numpoints.

plots:-implicitplot(
     eval(
          [[x^3+y^3-3*a*x*y, x+y+a], x= -3*a..3*a, y= -3*a..3*a],
          a= 1
     )[],
     gridrefine= 2, linestyle= [solid, dash], color= [red, blue]
);

Substitute z= exp(I*theta). Thus dz = I*exp(I*theta)*dtheta. 
z = 1 => theta(2*Pi); z = -1 => theta = Pi, and these limits trace the bottom half of the circle.

int(I*exp(I*t)/sqrt(exp(I*t)), t= 2*Pi..Pi);

The Euler totient is emphatically NOT the number of prime factors; indeed, it's nearly the opposite of that. The totient of n > 0 is the count of k, 0 < k <= n, such that k and n share NO prime factors. The symbol phi is usually used for this function. It's obvious from this definition that for prime p, phi(p) = p-1.

If n = p*q with p and q distinct primes, then phi(n) = (p-1)*(q-1). This fact is a very important aspect of RSA cryptograhy. Usually the only ways to compute the totient of a number all require the number's prime factorization. Since that is currently generally computationally infeasible when p and q are both very large, computing phi(n) is likewise infeasible. However, it hasn't been proven that there's no feasible algorithm. Quantum computation may be able to crack it if and when quantum computers get large enough.

Why don't you want to use Heaviside? A similar alternative is 

f1:= x-> piecewise(x >= 0, 1, 0);

If you want a more-primitive procedure based on if statements, then you need to delay the evaluation of x >= 0 in those cases where it can't be decided (such as when x has no value). Christian showed one way of doing that; here's another:

f1:= proc(x)
local r:= PiecewiseTools:-Is(x >= 0);
   `if`(r::truefalse, `if`(r, 1, 0), 'procname'(x))
end proc:

There's nothing wrong with your Maple: It can do the integral fairly quickly if you tell it, essentially, that the expression under the square root is nonnegative over the interval of integration. The conditions d/L__1 >= 2, d > 0 are sufficient. Use an assuming clause, like this:

simplify(...) assuming d/L__1 >= 2, d > 0;

where ... is exactly what you already have in the parentheses.

If these assumptions are not actually true for your problem, some deeper thought will be required.

It's nearly impossible to accurately compute the eigenvalues of a large matrix by finding the roots of the characteristic polynomial. Instead, use command LinearAlgebra:-Eigenvalues.

I show how to do that in this Post from a year ago: "Numerically solving BVPs that have many parameters" The BVP in that Post is extremely similar to yours; indeed, yours may be a direct simplififation of it. However, I wouldn't recommend that you attempt this if this is your first experience with Maple.

Like almost all operators in Maple, square brackets [] can be overloaded to have a special meaning for specific types of operands. Here is a simple implementation of what you're asking for:

Boolean_01:= module()
option package;
export `[]`:= proc(e::boolean, $) option overload; `if`(e,1,0) end proc;
end module:

with(Boolean_01);

#examples:
[3=7];
[3<7];

This is a minimalist implementation. There's no limit to the amount of subtlety that can be included.

1 /~ [8, 9, 9, 7, 9, 10, 5] mod 11

In general, to perforn an elementwise operation, append ~ to the operator. See help page ?elementwise.

The option scaling= constrained is the antithesis of being able to change the size of the z-axis. You should remove it.

Each seq command can only have one index variable. So, for example, if you wanted to generate all pairs of positive integers less than 4, you could use a double seq, like this:

seq(seq([i,j], i= 1..3), j= 1..3)

You could just as well iterate the j in the inside seq and i in the outside, which would generate the same pairs in a different order.

But you CAN'T do this:

seq([i,j], i= 1..3, j= 1..3)

simply because the syntax won't allow it, even though the concept is sensible. You attempted to do something like this, which was the reason for your error.

This is what I would do:

U:= eval(Array(1..2, 1..2, 1..2, 1..4, ()-> u[args]), Sol);

Now, U[2, 1, 1, 3], for example, would be the Array entry whose value corresponds to u[2, 1, 1, 3] in Sol. Unlike the other solutions based on either the assignment operator or assign command, this method does not cause any change to Sol or to the u's. Thus the u's can be used again as variables, perhaps to obtain another solution. This re-use of the u's would not affect the values now stored in the Array.

My method will work even if there are some index combinations for which there is no corresponding u in Sol.

My general policy is to avoid making any direct assignments to the fundamental variables of a problem (such as the independent and dependent variables); rather, I store solution values somewhere else, such as the Array, or just keep them in Sol for use on an as-needed basis. The command eval(expr, Sol) can be used to apply the values to any expression expr containing any of the u's.

In my code that you refer to, there is a line

Jxy:= ....

The purpose of that line is to list (in jet notation) the derivatives that you want to consider. In that case, I wrote some fancy code to extract all second derivatives of u and v with respect to x and/or y, but I could've just as well have written simply

Jxy:= {u[x,x], u[x,y], u[y,y], v[x,x], v[x,y], v[y,y]}.

Likewise, you can put whatever you want in the set. This is the only change that's needed.

Note that if the equation is nonlinear with respect to the desired derivatives, then a concept such as "the coefficient of u[x,x]" is not very meaningful. Nonlinearity of the overall equation is irrelevant; it's the linearity with respect to the desired derivatives that matters.

The problem is not recursive functions but a recursive assignment statement. Your first assignment statement inside the loop has x[n+1] on both sides of the assignment operator.