AbsoluteTiming[
 Parallelize[
  f[n_] := n (3 n - 1)/2;
  g[{i_, j_}] := IntegerQ[(1 + Sqrt[1 + 24 (f[i] + f[j])])/6] &&
    IntegerQ[(1 + Sqrt[1 + 24 (-f[i] + f[j])])/6];
  Select[Flatten[Table[{i, j}, {i, 1, 2000}, {j, i + 1, 3000}], 1],
   g]]]