Here is a way to recover the real solutions:
Clear[x1, x2, x3, x4, y1, y2, y3, y4, x, y, a, b, n];
{{x1[b_, n_], y1[b_, n_]}, {x2[b_, n_], y2[b_, n_]}, {x3[b_, n_],
y3[b_, n_]}, {x4[b_, n_],
y4[b_, n_]}} = {x, y} /.
Solve[{x^2 - y == a, y^2 - x == b}, {x, y}] /.
Solve[16^a + 16^b == 16, a][[1]] /. C[1] -> n;
ParametricPlot[{Chop@{x3[b, 0], y3[b, 0]},
Chop@{x4[b, 0], y4[b, 0]}}, {b, -2, 2}]
The result can be checked against this:
ContourPlot[16^(x^2 - y) + 16^(y^2 - x) == 16, {x, -1, 2}, {y, -1, 2}]