The direction field (by arrows) will only be displayed for one first-order equation or for an autonomous system of two first-order equations.

From the help for DEtools:-DEplot command: "A system of two first order differential equations produces a direction field plot, provided the system is determined to be autonomous. In addition, a single first order differential equation produces a direction field (as it can always be mapped to a system of two first order autonomous differential equations). A system is determined to be autonomous when all terms and factors, other than the differential, are free of the independent variable. For more information, see DEtools[autonomous]. For systems not meeting these criteria, no direction field is produced (only solution curves are possible in such instances). There can be only one independent variable."
 

Specifying a complex function in  n  complex variables is equivalent to specifying two real functions, each in  2*n  real variables.

Example for n=2 :

restart;
f:=(z1,z2)->(z1+z2)^2;
z1:=x1+I*y1: z2:=x2+I*y2:
u:=unapply(evalc(Re(f(z1,z2))),(x1,y1,x2,y2));
v:=unapply(evalc(Im(f(z1,z2))),(x1,y1,x2,y2));

(u+I*v)(1,2,3,4); # Example of a calculation

\

This expression is a real number (if  z  and  b  are real numbers):

A:=(ln (-z * I + b) - ln (z * I + b)) * I;
simplify(evalc(A)) assuming z>0,b>0 ;
simplify(evalc(A)) assuming z>0,b<0 ;

F:=beta[1]*exp(x*alpha[1]+y*beta[1]-z*sqrt(-alpha[1]^2-beta[1]^2))/(1+exp(x*alpha[1]+y*beta[1]-z*sqrt(-alpha[1]^2-beta[1]^2)));
evalc(F); # F in algebraic form (assuming all the parameters are reals)
evalc(Re(F));
evalc(Im(F));

I interpreted OP's question as how to find tetrahedra with integer edges and volume, and such that all edges are different. The code below finds all substantially different solutions, i.e. such tetrahedra, which cannot be obtained from one another by isometry, and whose edges are coprime, If the edge lengths <= 20, we obtain 68 solutions under these conditions. Also, all solutions are sorted in ascending order of volume (the first number in each list).

Edit. The code has been fixed. The necessary condition  X[1]=Y[1]  added to the  ListTools:-Categorize command.

Download itetr2.mw

Alternatively, you can do the same without matrices. Also, I multiplied each frame 5 times to make the animation slower.

T:= seq(k, k= 1..9):  X:= seq(k^2, k= 1..9):
plots:-display(seq(plot([seq([T[i],X[i]],i=1..k)])$5,k=1..9),insequence);

Use  mtaylor  command for this. 

Example:

restart;
i2 := (x,y) -> -(1/2)*I*(exp(I*x)*(sin(x)/x)-exp(I*y)*(sin(y)/y))/(x-y);
i3_r := -(1/2)*I*(i2(y,z)-i2(y,x))/(z-x);
mtaylor(i3_r, {x,y,z}, 6);

Because these 3 points are located in the same plane (the plane  yOz ), then you can use the  CurveFitting:-LeastSquares  command. We use a 1st degree polynomial  z=a*y+b . Increasing the degree gives nothing, because you have only 3 points, and two of them lie on the same vertical. A polynomial of any degree  z=P(y)  found by the least squares method will pass through point  C  and the middle of the segment  AB :

Download 3_points.mw

1. I made some corrections to your worksheet that Thomas Richard mentioned already.

2. I am not quite sure what exactly you want to plot. For each value of  y , equations  a  and  b  define some lines on the plane. Below on the first graph these lines are plotted for  y = 1 . On the second graph, using the  Explore  command, you can trace how these graphs change with the change in the variable  y .

Download eq_(1)_new.mw

Should be the  plots:-contourplot  command instead of the  contourplot one,  because it is from a pluggable package plots.

Remove the extra letter from the command name. Should be  plot3d   instead of  plotd3d .

The plot clearly shows: one maximum, two minimums and three saddles

Download zadelpunten_procedure_new.mw

Use the  seq  command for this.
Example:

restart;
L:=[0, Pi/6, Pi/4, Pi/3, Pi/2];
[seq([p,sin(p)], p=L)];

restart;
P:=x^4 - x^2 + x - 1:
A:=PolynomialTools:-Split(P, x);
convert(A, radical);
simplify(expand(%));  # A check

1. For  isolve a system should be specified as a set (not a list).
2. You need to specify exactly what you mean by the smallest positive solution, because the solution will not be one number, but an ordered set of numbers (a list). For example, which is less: [2, 4] or [1, 5] ?

In the solution below, the solutions are sorted by the minimum sum of the squares of all coordinates.

restart;
sys:= {a[1] = 2 * a[5], a[1] = a[4], 4 * a[1] + 2 * a[3] = a[6] + 2 * a[7], 2 * a[2] + 2 * a[3] = 2 * a[6], a[2] = a[4] + a[5], 2 * a[2] -a[1] = 2 * a[4]};
indets(sys);
n:=nops(%);
isolve(sys);
assign(%);
solve({seq(a[i]>0, i=1..n)});
[seq(seq([seq(a[i],i=1..n)],_Z2=floor(5*_Z1/2)+1..10),_Z1=1..10)];
sort(%, key = (x->add(x[i]^2,i=1..n)));