Such triangles exist even in the plane. (AE is the angle bisector in  ΔABC.)

restart;

IntegerTriangle := proc(iter)
local p,q,r,s, A:=[0,0], B,C,E,AB,AC,BC,AE,  den,
      nor:=V -> sqrt(V[1]^2+V[2]^2):
for p to iter do
for q to p-1  do   if igcd(p,q)>1 or irem(p,2)=irem(q,2) then next fi;
for r from p to iter do
for s to r-1  do   if igcd(r,s)>1 or irem(p,2)=irem(q,2) then next fi;
  B:=[p^2-q^2, 2*p*q]: AB:=p^2+q^2:
  C:=[r^2-s^2, 2*r*s]: AC:=r^2+s^2:
  BC:=nor(B-C):
  if not type(BC, rational) then next fi:
  E:=(AC*~B+AB*~C)/~(AC+AB):
  AE:=nor(E):
  if not type(AE, rational) then next fi:
  if (AB=AC)or(AB=BC)or(AC=BC) then next fi:
  if (AB+AC+BC)<=2*max(AB,AC,BC) then next fi:
  if (AB^2=AC^2+BC^2)or(AC^2=AB^2+BC^2)or(BC^2=AB^2+AC^2) then next fi:
  den:=ilcm(denom~([E[],AE])[]): B:=B*den: C:=C*den: E:=E*den;
  print(ABCE=[A,B,C,E]);
od od od od:
print("OK")
end proc:
IntegerTriangle(100);

 ABCE = [[0, 0], [4797, 1804], [4797, 280604], [4797, 6804]]
 ABCE = [[0, 0], [697, 14280], [18377, 14280], [7425, 14280]]
ABCE = [[0, 0], [1989, 116348], [55029, 45628], [34845, 72540]]
ABCE = [[0, 0], [7605, 444860], [226395, 153140], [143136, 264152] ]
"OK"

restart;
eq := x^2 + y^2 + z^2 - 2*x + 8*y - 6*z - 30 = 0:
Student[Precalculus][CompleteSquare](eq);

                            (z-3)^2 + (y+4)^2 + (x-1)^2 - 56  =  0

You must unevaluate the argument of latex using right single quotes 'expr'.

latex('diff(phi(t), t) = (`𝒟`[1, 4](phi(t)) + `𝒟`[4, 1](phi(t)))*diff(phi(t), t, x) + (`𝒟`[2, 4](phi(t)) + `𝒟`[4, 2](phi(t)))*diff(phi(t), t, y) + (`𝒟`[3, 4](phi(t)) + `𝒟`[4, 3](phi(t)))*diff(phi(t), t, z) + (`𝒟`[1, 2](phi(t)) + `𝒟`[2, 1](phi(t)))*diff(phi(t), x, y) + (`𝒟`[1, 3](phi(t)) + `𝒟`[3, 1](phi(t)))*diff(phi(t), x, z) + (`𝒟`[2, 3](phi(t)) + `𝒟`[3, 2](phi(t)))*diff(phi(t), y, z) + `𝒟`[4, 4](phi(t))*diff(phi(t), t, t) + `𝒟`[1, 1](phi(t))*diff(phi(t), x, x) + `𝒟`[2, 2](phi(t))*diff(phi(t), y, y) + `𝒟`[3, 3](phi(t))*diff(phi(t), z, z) + (diff(phi(t), t) + diff(phi(t), x) + diff(phi(t), y) + diff(phi(t), z))*(D(`𝒟`[1, j])(phi(t))*diff(phi(t), x) + D(`𝒟`[2, j])(phi(t))*diff(phi(t), y) + D(`𝒟`[3, j])(phi(t))*diff(phi(t), z) + D(`𝒟`[4, j])(phi(t))*diff(phi(t), t))');

Otherwise your expression is evaluated, and for example diff(phi(t), t, x) evaluates to 0.

A way to obtain the plot file is the following:

In the worksheet "f.mw" subject to RunWorksheet, don't use plotsetup but use:

...
p := plot(...):
...
return p;

(so, return the PLOT structure).
In the calling worksheet use:

ret := DocumentTools :-RunWorksheet("f.mw"):
plottools:-exportplot( "d:/temp/myplot.png", ret);

General solution: p, q, r, n arbitrary integers (n>0) and s = k^3 - p*A - q*B - r*C, where k is an arbitrary integer.
(Not very interesting).

Download Primfield-vv.mw

Just replace
S := sort([sqrt(x2), sqrt(y2), sqrt(z2)]);
with

S := sort([sqrt(x2), sqrt(y2), sqrt(z2)], key=evalf);

This is necessary because we need a numeric rather than symbolic ordering (see ?sort).

Q5:=n -> irem(add(convert(n,base,10)),5):
for n to 100000 do
if (Q5(n)+Q5(n+1)=0) then print(nmin=n); break fi od:

                            nmin = 49999

Download birth-vv.mw

Download sol-vv.mw

dsolve declares alias(c__1 = _C1) ;

restart;
dsolve(diff(y(x),x)=1);
alias();
_C1;

            y(x) = x + c__1
            c__1
            c__1

For container structures (such as tables, rtables) use e.g. EqualEntries  (or LinearAlgebra:-Equal(...), or ...)

EqualEntries(Matrix([[1,2],[1,2]]), <1,2;1,2>);

        true

The equation has a unique real solution, so, "Diophantine" is not essential.

restart;
f:=419*x^2 + 116*x*y - 426*x*z + 78*y^2 - 142*y*z + 133*z^2 - 1604*x - 682*y + 1086*z + 2306:
CS:=Student:-Precalculus:-CompleteSquare:
CS(f,z):
f1:=select(has,[op](%), z)[]:
CS(f-f1,y):
f2:=select(has,[op](%), y)[]:
f3:=CS(expand(f-f1-f2),x):
expand(f-f1-f2-f3);
solve([f1,f2,f3]);

          0
          {x = 7, y = 11, z = 13}

By Dirichlet's theorem on arithmetic progressions - Wikipedia there are infinitely many primes of the form 24*n + 1.
p,q can be any pair of these.

P.S. We could avoid Dirichlet' theorem noting that at least one of the sets {24*n+k| n in N} for k=0,1,...,23   must contain infinitely many primes.

You want the limit of V for beta --> oo. Use:

limit(V, beta = infinity)

It results 6.5.

If you want to check this numerically for large values for beta (as you did), increase first Digits e.g.

Digits := 30;

P.S. Note that usually, for symbolic computations like limit, it is recommended to use exact values or rational values for parameters.