I. Third Powers
As background, the complete solution to $x_1^3+x_2^3 = x_3^3+x_4^3$ was given by Euler and Binet using deg-$4$ polynomials. However, Elkies also found a complete solution using deg-$3$ polynomials in $3$ variables. Define.
$$f(a,b,c) = 9a^3 + 9a^2b + 3a^2c + 3a b^2 - 6a b c + 3a c^2 + 3b^3 + 3b^2c + b c^2 + c^3$$
then,
$$f(a,b,c)^3+f(-a,b,-c)^3=f(a,b,-c)^3+f(-a,b,c)^3$$
II. Fourth Powers
In contrast, no complete polynomial solution to $x_1^4+x_2^4 = x_3^4+x_4^4$ is known (or maybe even possible). However, analogous to Elkies' solution, there are at least five of form,
$$f(a, b)^4 + f(b, -a)^4 = f(a, -b)^4 + f(b, a)^4$$
where $f(a,b)$ is a polynomial of degree $k = 6n+1 = 7,13,19,25,31$ given in the addendum below. Whether there are more remains to be seen. (Deg-$25$ was found by yours truly.)
III. Elliptic curve
Using Euler's approach, let,
$$(p+q)^4+(r-s)^4=(p-q)^4+(r+s)^4$$
and define $p = (a^3 - b),\, q = a y,\, r = b (a^3 - b),\, s = y.\,$ Then the equation above transforms to the simple form,
$$(a^3 - b) (b^3 - a) = y^2$$
This is birationally equivalent to an elliptic curve. (Update: This curve has been analyzed by Deyi Chen in this MO answer.) Known rational points of small height are,
$$a=n,\quad b =n\,\frac{(4 + n^2 + 10n^4 + n^6)}{(1 + 10n^2 + n^4 + 4n^6)}\quad$$
$$a=n^3,\quad b =n\,\frac{(1 - 2n^2 + n^4 + n^6)}{(1 + n^2 - 2n^4 + n^6)}\quad$$
$$\quad a=n,\quad b =n\,\frac{(9 - 44n^2 + 190n^4 + 100n^6 + n^8)}{(1 + 100n^2 + 190n^4 - 44n^6 + 9n^8)}$$
$$\quad a=n^3,\quad b =n\,\frac{(1 - 4n^2 + 10n^6 + 8n^{10} + n^{12})}{(1 + 8n^2 + 10n^6 - 4n^{10} + n^{12})}$$
Note the symmetry of $b$'s numerator and denominator. The first two points $b$ yield the same deg-$7$ solution (after removing common factors) while the last two $b$ yield the deg-$13$ and deg-$19$ solutions. For the deg-$25$, we have $a=n,$ and,
$\qquad\qquad b =\frac{n(16 - 543n^2 + 4632n^4 + 15100n^6 + 10632n^8 + 22758n^{10} + 6568n^{12} + 5820n^{14} + 552n^{16} + n^{18})}{(1 + 552n^2 + 5820n^4 + 6568n^6 + 22758n^8 + 10632n^{10} + 15100n^{12} + 4632n^{14} - 543n^{16} + 16n^{18})}$
found by yours truly (though I may have missed a point of smaller height). The deg-$31$ solution has a rational point $b$ that is of higher height.
IV. Question
Q: Given some fixed $a$ (like $a=n$) and the elliptic curve $(a^3 - b)(b^3 - a) = y^2,$ can you find a rational point $b = P(n)/Q(n)$ that is of height less than deg-$24$?
Note: Points that yield an identity of form $f(a, b)^4 + f(b, -a)^4 = f(a, -b)^4 + f(b, a)^4$ are preferred. An example of a point that does not do so was found by Lander (and which leads to a second deg-$19$ identity).
Addendum:
Below are the five known polynomials $f(a,b)$ of deg $k = 6n+1 = 7,13,19,25,31.$ (For ease of copy-paste, they are purposely not in Latex format for those who want to test them.)
The deg-$7$, or versions thereof, have been independently found by many including Euler, Gerardin, Swinnerton-Dyer (at the tender age of 16), and others.
Deg 7
f(a, b) = a^7 + a^5 b^2 - 2 a^3 b^4 - 3 a^2 b^5 + a b^6
Deg 13
f(a, b) = a^13 + a^12 b - a^11 b^2 + 5 a^10 b^3 - 6 a^9 b^4 - 12 a^8 b^5 + 4 a^7 b^6 + 7 a^6 b^7 + 3 a^5 b^8 - 3 a^4 b^9 - 4 a^3 b^10 + 2 a^2 b^11 + a b^12 + b^13
Deg 19
f(a, b) = a^19 + 6 a^17 b^2 - 18 a^15 b^4 + 6 a^14 b^5 - 5 a^13 b^6 + 12 a^12 b^7 - 12 a^11 b^8 + 36 a^10 b^9 - 24 a^9 b^10 - 12 a^8 b^11 + 19 a^7 b^12 + 36 a^6 b^13 + 6 a^5 b^14 + 12 a^4 b^15 - 6 a^3 b^16 + 6 a^2 b^17 + a b^18
Deg 25
f(a, b) = a^25 + 7 a^23 b^2 - 2 a^21 b^4 - 3 a^20 b^5 - 25 a^19 b^6 - 63 a^18 b^7 + 43 a^17 b^8 + 36 a^16 b^9 - 134 a^15 b^10 + 213 a^14 b^11 + 179 a^13 b^12 - 333 a^12 b^13 - 128 a^11 b^14 + 207 a^10 b^15 + 136 a^9 b^16 - 57 a^8 b^17 - 97 a^7 b^18 + 9 a^6 b^19 + 4 a^5 b^20 - 9 a^4 b^21 + 10 a^3 b^22 - 3 a^2 b^23 + a b^24
Deg 31
f(a, b) = a^31 - a^30 b + 11 a^29 b^2 + a^28 b^3 + 42 a^27 b^4 + 24 a^26 b^5 - 19 a^25 b^6 - 32 a^24 b^7 - 154 a^23 b^8 - 254 a^22 b^9 + 266 a^21 b^10 + 718 a^20 b^11 + 126 a^19 b^12 - 303 a^18 b^13 - 478 a^17 b^14 - 830 a^16 b^15 + 770 a^15 b^16 + 916 a^14 b^17 - 738 a^13 b^18 + 21 a^12 b^19 + 350 a^11 b^20 - 434 a^10 b^21 + 50 a^9 b^22 + 142 a^8 b^23 - 91 a^7 b^24 + 76 a^6 b^25 + 15 a^5 b^26 - 3 a^4 b^27 + 8 a^3 b^28 - 8 a^2 b^29 + a b^30 - b^31
P.S. A similar question was asked by emacs in this post.
I am trying to go structurally, details will be in there so that the yoga is transparent at lowest level. (Instead, a computer one-liner does often the job, and this is not only frustrating for the reader, but also for me when years later i do not recall the human argument and regret that typing in hurry.) Let us start.
Consider the curve defined by $$ v^2 = (u^3-A)(A^3 -u)\ , $$ first as a(n affine) curve in the variables $(u,v)$ over the field $F=\Bbb Q(A)$. So $A$ is felt like a fixed parameter. We bring it in the short Weierstrass form using the substitution $\displaystyle u=A^3-\frac 1{x_1}$, and $\displaystyle v=\frac 1{x_1^2}y_1$, this leads to $$ \begin{aligned} y_1^2 &= (A^3x_1 -1)^3 - Ax_1^3\ ,\\ y_1^2 &= \underbrace{(A^9-A)}_{:=c}x_1^3 - 3A^6x_1^2 +3A^3x_1 - 1\ ,\\ (cy_1)^2 &= (cx_1)^3 -3A^6(cx_1)^2 +3A^3c(cx_1)-c^2\ . \end{aligned} $$ The coefficient in $x_1^3$ on the R.H.S was $c=(A^9-A)$ in the second line. So we multiplied with $c^2$ on both sides, then used $y_2=cy_1$, $x_2=cx_1$. Here is the transformed equation: $$ y_2^2 = x_2^3 - 3A^6x_2^2 + 3A^3cx_2-c^2\ ,\qquad c=A^9-A\ . $$ The first two terms on the R.H.S. suggest the substitution $x=x_2-A^6$. We take $y=y_2$ and obtain a pretty simple form: $$ \bbox[yellow]{\qquad (E=E(A))\ :\qquad\qquad y^2 = x^3 -3A^4x-A^2(A^8+1)\ . \qquad } $$ When we pass from $(u,v)$ to $(x,y)$ the connecting formulas make $u$ depend only on $x$, and conversely. Explicitly: $$ \begin{aligned} x &= x_2-A^6=cx_1-A^6 =\frac c{A^3-u} - A^6\ ,\\ u &= A^3-\frac 1{x_1}=A^3-\frac c{x_2}=A^3-\frac c{x+A^6}\ . \end{aligned} $$ Do we have a quick guess of a point on $E(A)$? We indeed have such a guess if we specialize (base change, composition with $\Bbb Q[A]\to \Bbb Q[n]$, $A\to n^3$) to the case $A=n^3$. Then the factor $(u^3-A)=(u^3-n^3)$ has an easy point, $u=n$, and we obtain a point $H$ with $H_x=c/(A^3-n)-A^6=(A^9-A)/(A^3-n)-A^6=(n^{27}-n^3)/(n^9-n)-n^{18} =n^2(n^8+1)$. The corresponding $y$-value is $H_y=0$ (since the vanishing of $(u^3-A)$ makes $v$ vanish, which propagates to a vanishing $y$-value). Of course, the obtained point $$ H = H(n^3) = (H_x, H_y) = (n^2(n^8+1),\ 0) $$ on $E(n^3)$ is a two-torsion point. We need "more points", other (kind of) points.
What about the general case? It turns out that there is also a simple point on the $E(A)$ curve, namely $$ \bbox[yellow]{ \qquad G =G(A) = (G_x,G_y)=(A^4+A^2+1,\ A^6+A^4+A^2+1)\ . \qquad} $$ In general, when we specialize $A$ to a rational number, there is "only this" point (with its multiples) in the list of $\Bbb Q$-rational points. So we cannot expect a bigger rank for $E(A)$ over $\Bbb Q(A)$.
With the above arguments, it is easy to see why the cases
lead to solutions for the initial problem. In the context of the OP the following questions arise:
Question A: Which are explicitly the solutions? For $A=n$ and for $A=n^3$.
Question B: What about "specializations" of $A$ to rational numbers? In which cases do we get a "high rank"? Can we guess a formula for a generic family $A=\mathcal A(n)$ so that new parametric points appear?
Question C: What about "other generic specializations" $A=\mathcal A(n)$ possibly coming from B?
Answer to Question A, when parameter $A$ is $A=n$.
We start with $G$. Let $k$ be $1,2,3,4,\dots$ and compute $kG=(x_k,y_k)$. We need only $x_k$. Compute the corresponding $b_k =u_k:=A^3 - c/(x_k+A^6)$. It turns out that $b_k=b_k(n)$ is a function of the shape $b_k(n)=n\cdot\frac{g_k(n,1)}{g_k(1,n)}$.
Such a $b_k$ determines $v_k$ from $$ v_k^2 =(a^3-b_k)(b_k^3-a) \ . $$ We use it to compute $p,q,r,s$ as in the OP. We are interested in $p-q$, which is a polynomial in $n$. We compute it in $n=(N+1)/(N-1)$, and after clearing the denominator, a power of $(N-1)$, we obtain a polynomial $f_k=f_k(N)$. The postponed code gives the following tables for the $b$-coefficients, and polynomials $f$:
and respectively
Used code (with True, False as shown to get the $f$-polynomials, respectively with False, True instead to get the $g$-polynomials):
... to be continued ...
(There is a maximal amout of characters that can be used in an answer, and the tables take them all...)