This can be done in many ways. In such cases, I usually write the simplest procedures using standard formulas:

restart;
dist:=(A,B)->sqrt((A[1]-B[1])^2+(A[2]-B[2])^2):
midp:=(A,B)->(A+B)/2:

# Example
dmf1:=[-3, 1]:  dmf2:=[3, 2]:
dist(dmf1,dmf2);
midp(dmf1,dmf2);

If we introduce coordinates  C(0,0), D(486,0), A(x1,y1), B(x2,y2) , we get an example of an underdetermined problem, since we have 4 unknowns, but only 3 conditions (3 known distances) to find them. Such systems usually have an infinite number of solutions. To find a specific solution, it is necessary to assign a certain value to one unknown. I considered 2 cases: x2=120  and  x2=140 . To find the area, I used the Area procedure, a detailed review of it is available at the link  https://mapleprimes.com/posts/145922-Perimeter-Area-And-Visualization-Of-A-Plane-Figure-

restart;
local D:

Area := proc(L) 
local i, var, e, e1, e2, P; 
for i to nops(L) do 
if type(L[i], listlist(algebraic)) then 
P[i] := (1/2)*add(L[i, j, 1]*L[i, j+1, 2]-L[i, j, 2]*L[i, j+1, 1], j = 1 .. nops(L[i])-1) else 
var := lhs(L[i, 2]); 
if type(L[i, 1], algebraic) then e := L[i, 1]; 
if nops(L[i]) = 3 then P[i] := (1/2)*(int(e^2, L[i, 2])) else 
if var = y then P[i] := (1/2)*simplify(int(e-var*(diff(e, var)), L[i, 2])) else 
P[i] := (1/2)*simplify(int(var*(diff(e, var))-e, L[i, 2])) end if end if else e1 := L[i, 1, 1]; e2 := L[i, 1, 2]; 
P[i] := (1/2)*simplify(int(e1*(diff(e2, var))-e2*(diff(e1, var)), L[i, 2])) end if end if end do; 
add(P[i], i = 1 .. nops(L)); 
end proc:

C:=[0,0]: D:=[486,0]: A:=[x1,y1]: B:=[x2,y2]:
ds:=(M,N)->(M[1]-N[1])^2+(M[2]-N[2])^2:
eqs:=ds(C,B)=(162.8)^2, ds(B,A)=(71.8)^2, ds(A,D)=325^2;
sol1:=fsolve({eqs,x2=120}, {x1=0..infinity,y1=60..infinity,y2=0..infinity,x2=0..infinity});
sol2:=fsolve({eqs,x2=140}, {x1=0..infinity,y1=60..infinity,y2=60..infinity,x2=0..infinity});
L:=[[C,D,A,B,C]]:

Area(eval(L, sol1)); # First solution
Area(eval(L, sol2));  # Second solution

See https://math.stackexchange.com/questions/513071/tetrahedron-inequality?newreg=b9e402d644c3454d92b79b0dd7d32a25

Below in the code, for brevity, I limited myself to the first 10 prime numbers:

restart;
with(combinat):
primes := [seq(ithprime(i), i = 1 .. 10)];
S:=choose(primes, 6):
F:=(x,y,z,u,v,w)->`if`(LinearAlgebra:-Determinant(<0,1,1,1,1; 1,0,x^2,y^2,z^2; 1,x^2,0,u^2,v^2; 1,y^2,u^2,0,w^2; 1,z^2,v^2,w^2,0>)>0,true,false): k:=0:
for s in S do
for t in [seq(seq([s1[],s2[]], s2=permute(convert(s,set) minus convert(s1,set))), s1=choose(s, 3))] do
x,y,z,u,v,w:=t[];
if F(t[]) and min(x+y+z,x+v+u,w+y+u,w+v+z)>max(x+w,y+v,z+u) then k:=k+1; L[k]:=t; break fi;
od: od:
L:=convert(L, list); 
nops(L);

We received the list  L  of 22 solutions.

I've seen many cases where this command fails to work if the system contains inequalities. The workaround is obvious:

restart;
L:=[isolve(a*b=4)];
select(s->`and`(op(rhs~(s)>=~1)), L)[];

PS. For such cases, the  IntegerPoints2   procedure from the post  https://mapleprimes.com/posts/202542-Integer-Points-In-Curved-Regions is useful:

IntegerPoints2({a*b=4,a>=1,b>=1},[a,b]);

@nm  I see that you are more interested in plotting integral curves (stream lines) than just arrows. Below I show how to plot 10 curves for your example. It is not very difficult to plot arrows along each line as well, or just arrows, using the  plots:-arrow  command. The picture clearly shows that the origin is a stable focus.

restart;
with(plots):
sys:=[diff(x(t),t) = -3*x(t)-4*y(t), diff(y(t),t) = 2*x(t)+y(t)]:
S1:=seq(dsolve({sys[],x(0)=4,y(0)=i}), i=0..4):
S2:=seq(dsolve({sys[],x(0)=-4,y(0)=i}), i=-4..0):
P:=plot([seq(eval([x(t),y(t),t=0..2*Pi],s), s=[S1,S2])], x=-4.5..4.5, y=-5..5, color=red, thickness=2):
V1:=seq(seq(arrow(eval(eval([x(t),y(t)],S1[k]),t=T),eval(eval([diff(x(t),t),diff(y(t),t)],S1[k]),t=T),length=0.9,head_width=0.15,head_length=0.2,color=red,border=false), T=[0.05,0.15,0.3,0.6,1.2]), k=1..5):
V2:=seq(seq(arrow(eval(eval([x(t),y(t)],S2[k]),t=T),eval(eval([diff(x(t),t),diff(y(t),t)],S2[k]),t=T),length=0.9, head_width=0.15,head_length=0.2,color=red,border=false), T=[0.05,0.15,0.3,0.6,1.2]), k=1..5):
display(V1,V2,P, size=[800,800], labels=[``,``], font=[times,16]);

The answer has been significantly edited.

stream_lines.mw

The  legend  option is not supported when drawing arrows. Instead, you can make text labels near each vector. In my opinion, this looks even better than using legends:
 

Download legend_question_new.mw

The equation that results here resembles Pell's equation, and apparently has infinitely many solutions. The  isolve  command can't handle this equation (returns  NULL) . In double for-loop all solutions for  n <=10000  are found :

restart;
sum(i,i=1..k-1)=sum(i,i=k+1..n);
eq:=expand((lhs-rhs)(%)*2);
isolve(eq=0);
for n from 3 to 10000 do
for k from 2 to n-1 do
if eq=0 then print('k'=k,'n'=n) fi;
od: od: 

Here is my attempt at a solution. The solution is obtained independently of vv's one, although in essence both solutions coincide. The birds meet in the orthocenter of the tetrahedron, this point is equidistant from its vertices and is located at the height of the tetrahedron, dropped from the vertex at the distance of 1/4 of the height  from the base. I didn't find the "last" option in my version Maple 2018.2, probably it appeared only in the latest versions. The animation takes place in real time, approximately 8 seconds (birds fly at a speed of 0.1 m/sec, the side of the tetrahedron is 1 meter):

restart;
d:=v->sqrt(v[1]^2+v[2]^2+v[3]^2):
A:=[1/2,-sqrt(3)/6,0]: B:=[0,sqrt(3)/3,0]: C:=[-1/2,-sqrt(3)/6,0]: E:=[0,0,sqrt(6)/3]: 
fA:=t-> [xA(t),yA(t),zA(t)]: fB:=t-> [xB(t),yB(t),zB(t)]: fC:=t-> [xC(t),yC(t),zC(t)]: fE:=t-> [xE(t),yE(t),zE(t)]:
eq:={seq(diff(fA(t),t)[i]=0.1*(fB(t)-fA(t))[i]/d(fB(t)-fA(t)),i=1..3),seq(diff(fB(t),t)[i]=0.1*(fC(t)-fB(t))[i]/d(fC(t)-fB(t)),i=1..3),seq(diff(fC(t),t)[i]=0.1*(fE(t)-fC(t))[i]/d(fE(t)-fC(t)),i=1..3),seq(diff(fE(t),t)[i]=0.1*(fA(t)-fE(t))[i]/d(fA(t)-fE(t)),i=1..3),seq(fA(0)[i]=A[i],i=1..3),seq(fB(0)[i]=B[i],i=1..3),seq(fC(0)[i]=C[i],i=1..3),seq(fE(0)[i]=E[i],i=1..3)}:
sol:=dsolve(eq, numeric):

with(plots): with(plottools):
Tetr:=display(curve([A,B,C,A,E]),curve([B,E,C]), color=black):
P:=odeplot(sol,[[xA(t),yA(t),zA(t), color=red, thickness=3],[xB(t),yB(t),zB(t), color=blue, thickness=3],[xC(t),yC(t),zC(t), color=green, thickness=3],[xE(t),yE(t),zE(t), color=gold, thickness=3]],0..8.12, frames=80):
P1:=animate(s->display(line(eval([xA(t),yA(t),zA(t)],sol(s)),eval([xB(t),yB(t),zB(t)],sol(s)),linestyle=3,color=black)),[t],t=0..8.12,frames=80):
P2:=animate(s->display(line(eval([xC(t),yC(t),zC(t)],sol(s)),eval([xB(t),yB(t),zB(t)],sol(s)),linestyle=3,color=black)),[t],t=0..8.12,frames=80):
P3:=animate(s->display(line(eval([xC(t),yC(t),zC(t)],sol(s)),eval([xE(t),yE(t),zE(t)],sol(s)),linestyle=3,color=black)),[t],t=0..8.12,frames=80):
P4:=animate(s->display(line(eval([xA(t),yA(t),zA(t)],sol(s)),eval([xE(t),yE(t),zE(t)],sol(s)),linestyle=3,color=black)),[t],t=0..8.12,frames=80):
T:=textplot3d([[A[],"A",align=below],[B[],"B",align=below],[C[],"C",align=below],[E[],"E",align=above]], font=[times,17]):
display(Tetr,P,P1,P2,P3,P4,T, orientation=[45,75],axes=none);

Edit. 

The animation has been improved. I have added dash-lines connecting each bird to the bird it is flying towards. We can see that these lines are tangent to the corresponding trajectories.

birds.mw

Christmas1.mw
Let's introduce the coordinates. We'll take the length of each side as 1. We'll take E as the origin of the coordinates, we'll direct the x-axis along ED, and the y-axis vertically upwards. To calculate the lengths and angles, we'll use vector operations:

restart;
local D:
Digits:=20:
d:=v->sqrt(v.v):
E:=<0,0>: D:=<1,0>: C:=<1,1>: B:=<x,y>: A:=<cos(2*Pi/3),sin(2*Pi/3)>:
BA:=A-B: BC:=C-B:
sol:=simplify([solve({d(BA)=1,d(BC)=1}, explicit)])  assuming real;
evalf(%); # Coordinates of point B (only the first solution is suitable)
assign(sol[1]):
AE:=E-A: AB:=-BA: CB:=-BC: CD:=D-C:
a:=evalf(arccos(AE.AB)*180/Pi); b:=evalf((2*Pi-arccos(BA.BC))*180/Pi); c:=evalf(arccos(CB.CD)*180/Pi); # Angles in degrees
90+120+a+b+c;  # Check. The sum of the interior angles in any pentagon is 540 degrees.

Visualization:

with(plots): with(plottools):
L:=evalf(eval(convert~([E,A,B,C,D,E],list))):
P1:=curve(L, color=blue, thickness=3):
P2:=polygon(L, color="LightBlue"):
E,A,B,C,D:=L[1..5][]:
T:=textplot([[E[],"E",align=[left,below]],[A[],"A",align=[left,above]],[B[],"B",align=[right,below]],[C[],"C",align=right],[D[],"D",align=[right,below]]],font=[times,16]):
display(P1,P2,T, scaling=constrained, axes=none);

From your drawing we see that the vector F_RC1_ST is the opposite of the projection of F_RC1 on the z axis, and the vector F_RC1Y is equal to - F_RC1_ST - F_RC1 . In the figure below, all the vectors found come out from the origin. Figure out for yourself from the help how to draw them so that they form a triangle (to do this, you just need to change the bases of the two vectors from three ones).

>  >  EV_C1Y:=   >  (1) >

Download vector_question_new1.mw

Of course, the method @Carl Love showed is the most optimal, but it is always useful to find alternative ways. Imagine that after some time you need to solve a similar problem again, but you have completely forgotten these commands. In fact, the problem is easily solved iteratively by the simplest for-cycle for a list of any length:

restart;
lst1 := [a^3/10, -a^2/2, a, a^4/4]:
n:=nops(lst1): d[1]:=gcd(lst1[1..2][]):
for i from 3 to n do
d[i-1]:=gcd(d[i-2],lst1[i]);
od:
d[n-1];

If you need to do this multiple times, the best way to write a procedure for this:

restart;
P:=x->parse(convert(x,string)[1]):
P(asdf);

The imaginary unit in Maple is encoded as  I  rather than i . Maple performs the simplest arithmetic operations on complex numbers automatically, i.e. there is no need for any commands:

(3-4*I)*(2+I);

I don't know any other way to solve the problem than to do it using the tools of the plots and plottools packages:

restart; 
with(plots): with(plottools): 
P := plot(sin(x), x = -3 .. 3, colour = [blue], style = pointline, symbol = solidcircle, numpoints = 20): 
p1 := line([1.5, -1], [1.9, -1], colour = [blue]): 
p2 := pointplot([1.7, -1], symbol = solidcircle, color = blue, symbolsize = 10): 
p3 := curve([[1.4, -.9], [2.65, -.9], [2.65, -1.1], [1.4, -1.1], [1.4, -.9]]): 
p4 := textplot([2.3, -1, sin(x)]): 
plots:-display(P, p1, p2, p3, p4);

Download Plot1_new.mw

The problem can be easily reduced to solving a system of 2 equations with 2 unknowns. But for some reason, there is no correct answer among the suggested ones (possibly a typo):

restart;
d:=(X,Y)->sqrt((X[1]-Y[1])^2+(X[2]-Y[2])^2):
A:=[7,6]: B:=[3,4]: P:=[x,0]: Q:=[0,y]:
solve({d(P,A)=d(P,B),d(Q,A)=d(Q,B)});
assign(%);
d(P,Q); # The answer