quartals[m_, p_, q_, max_] :=
Module[{ans = {}, lhsD = <||>, lhs, v, u, w, rhs},
For[u = 1, u <= max, u++, lhs = m^p + u^q;
AssociateTo[lhsD, lhs -> Append[Lookup[lhsD, lhs, {}], u]]; ];
For[v = m + 1, v <= max, v++,
For[w = 1, w <= max, w++, rhs = v^p + w^q;
If[KeyExistsQ[lhsD, rhs],
Do[AppendTo[ans, {m, u, v, w}], {u, lhsD[rhs]}]; ]; ]; ];
ans = SortBy[ans, #[[2]] &];
Do[Print["Solution ", i, ": ", ans[[i]], " (", m, "^", p, " + ",
ans[[i, 2]], "^", q, " = ", ans[[i, 3]], "^", p, " + ",
ans[[i, 4]], "^", q, " = ", m^p + ans[[i, 2]]^q, ")"], {i,
Length[ans]}]; ans];
solns = quartals[2, 1, 3, 2000] (* solutions restricted to v<2000 *)
Grid[solns]
u1 = Map[#[[2]] &, solns] (*u,
A003057 *)
v1 = Map[#[[3]] &, solns] (*v,
A386215 *)
w1 = Map[#[[4]] &, solns] (*w,
A004736 *)