Try changing the initial conditions slightly by taking  y(0)=0.001  instead  y(0)=0. If  y(0)=0, then the vector of derivatives will be equal to  0  and the solution cannot be continued:

Download DEtools_example_new.mw

Here is the Maple solution to this problem. For the solution, 3 my procedures Perimeter, Area and Picture were used (for details see  https://www.mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure- ). For convenience, procedure codes are inserted into the document below. The idea behind the solution is to replace this shape with a polygon with a large number of vertices (2177 ones):

Download 3_proc.mw

@The function The example  b)  was solved incorrectly. You are trying, following Carl's instructions, to express the modulus and argument of a complex number  z  through its real and imaginary parts, but you are doing it wrong. It will be much easier to use direct commands  abs  and  argument  for  z  in conjunction with  evalc  command, which tells Maple that  x  and  y  are real numbers. The empty circle (in the pic below) means that the origin does not belong to the boundary of this region:

restart;
z:=x+I*y:
A:=plots:-inequal(evalc~({argument(2*I*z)>=Pi/2,argument(2*I*z)<=Pi,abs(4*I*z)>0,abs(4*I*z)<=1}), x=-0.35..0.35, y=-0.35..0.35, thickness=2, optionsfeasible=[color="LightBlue"]):
B:=plots:-pointplot([[0,0]], symbol=circle, symbolsize=18):
plots:-display(A,B, size=[500,500]);

When you write an arbitrary complex number  z  in algebraic form   z=x+I*y , Maple does not know that you are considering  x  and  y  as real numbers. Therefore, the commands  abs(z)  or  argument(z)  do not work directly. In this case, you can use the  evalc  command, which assumes that all symbols are real numbers. Also you can write  assuming real  instead. Then you do not need to write formulas for abs(z)  or  argument(z)  yourself, Maple will do it automatically. Below is the example  d)  solution:

restart;
z:=x+I*y:
abs(z);
evalc(abs(z));
plots:-inequal({abs(z)>=1,abs(z)<3,argument(z)<3*Pi/4,argument(z)>=-Pi/4}, x=-4..4, y=-4..4) assuming real;

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 .