It is actually quite straightforward to write down examples in one variable where this occurs. For example, the Diophantine equation $(x^2 - 2)(x^2 - 3)(x^2 - 6) = 0$ has this property: for any prime $p$, at least one of $2, 3, 6$ must be a quadratic residue, so there is a solution $\bmod p$, and by Hensel's lemma (which has to be applied slightly differently when $p = 2$) there is a solution $\bmod p^n$ for any $n$. We conclude by CRT. (Edit: As Fedor says, there are problems at $2$. We can correct this by using, for example, $(x^2 - 2)(x^2 - 17)(x^2 - 34)$.)
Hilbert wrote down a family of quartics with the same property. There are no (monic) cubics or quadratics with this property: if a monic polynomial $f(x) \in \mathbb{Z}[x]$ with $\deg f \le 3$ is irreducible over $\mathbb{Z}$ (which is equivalent to not having an integer solution), then by the Frobenius density theorem there are infinitely many primes $p$ such that $f(x)$ is irreducible $\bmod p$.
Oh, I think the answer is definitely yes!
Let $\{k \to x,y\}$ be any solution of $x^2 - (k^2+1)y^2 = k^2$, and let $K$ be the set of $k$ for which a solution has $0 < k < y-1$. In a paper recently
submitted to Glasnik Matematicki we call these solutions
exceptional solutions. Andrej's conjecture is that for any $k$ there is at most 1 exceptional solution.
One interesting result we obtain is that, if $k \in K$, then
$y < (2 - \sqrt{3})k$.
A particular feature of this Pell equation is its symmetry wrt $k$ and $y$. These are interchangeable, so for any solution $\{k \to x,y\}$ there is a corresponding solution $\{y \to x,k\}$.
It follows that if $k \neq y \pm{1}$, then either $k \in K$ or $y \in K$.
Now, for any $k \geq 2$, we have 3 particular solution classes $(x_n, y_n)$ with
$y_0 = \{0, k-1, -(k-1)\}$. For any $n > \{0, 0, 1\}$ we have $y_n > k-1$
and so $\{y_n \to x_n,k\}$ is exceptional, ie. $y_n \in K$.
We also need to consider $k=1$, for which there is just the one
class with $(x_0, y_0) = (1, 0)$ , $(x_1, y_1) = (3, 2)$ and so for
any $n > 1$ we have $y_n \in K$. For example $(x_2, y_2) = (17, 12)$ from
which we obtain exceptional solution $\{12 \to 17, 1\}$, and so $12 \in K$.
In our paper we call the corresponding set of exceptional solutions
"Type 1". But here let us simply define the set $K_1$ to be all of these
$y_n > k$ that we find from these 3 classes for any $k > 1$, and
from the one class for $k = 1$.
One property shared by all type 1 solutions, ie $\{k \to x, y\}$ with $k \in K_1$ is that either $(y^2 +1) | (x+y)$ or $(y^2 + 1) |(x-y)$.
Now, for any $k \in K_1$ we have a corresponding $\{k \to x,y\}$ for
which our Pell eqn has 2 additional solution classes, with fundamental solutions
$(x_0, y_0) = (x, \pm y)$. For any $n > \{0,1\}$ we then have $y_n > k-1$
and so $\{y_n \to x_n, k\}$ is exceptional, ie $y_n \in K$.
And of course we can apply the same process to any of these new $y_n$
ad infinitum, each $y_n$ seeding a forest of others. For example,
just considering $n = 1$ alone in each case, from $\{8 \to 18,\pm{2}\}, 8 \in K_1$
we obtain $\{546 \to 4402,8\}$ and $\{30 \to 242,8\}$, so $546, 30 \in K$,
and from $\{30 \to 242, \pm{8}\}$ we get $\{28928 \to 868322,242\}$ and
$\{112 \to 3362,30\}$, so $28928, 112 \in K$.
We call these "Type 2" solutions, so let's define $K_2$ to be all of
the $y_n$ found this way. These do not have the divisibility property that
was noted above for the $y_n \in K_1$.
In the paper we show that all
exceptional solutions can be enumerated recursively in this
fashion, ie. that $K = K_1 \cup K_2$. This is done by showing any $k \in K$ can
be traced back to a root in $K_1$.
The enumeration
algorithm is given below. Solution classes are referred to as $0, -1, +1$, the interpretation of which I hope is reasonably clear!
Proc Enum_K:
Enum_K1(1,0)
for k = 2 to $ \infty $
Enum_K1(k, 0)
Enum_K1(k, +1)
Enum_K1(k, -1)
Proc Enum_K1(k, class):
set $(x_0, y_0), (x_1, y_1)$ according to class
n1 = if (class = -1 or k = 1) then 2 else 1
for n = n1 to $\infty$
add $y_n$ to $K_1$
Enum_K2($y_n$, +1)
Enum_K2($y_n$, -1)
Proc Enum_K2(k, class):
set $(x_0, y_0), (x_1, y_1)$ according to class
n1 = if (class = -1) then 2 else 1
for n = n1 to $\infty$
add $y_n$ to $K_2$
Enum_K2($y_n$, +1)
Enum_K2($y_n$, -1)
To generate the solution sequences in any class, we note that each class has the same
recurrence relation:
$R = 2k^2 + 1$
$x_n = 2Rx_{n-1} - x_{n-2}$
$y_n = 2Ry_{n-1} - y_{n-2}$
but of course have different initial conditions:
$R = 2k^2 + 1, S = 2k, D = k^2 + 1$
$K_1, class 0: (x_0, y_0) = (k, 0), (x_1, y_1) = (kR, kS)$
$K_1, class +: (x_0, y_0) = (k^2-k+1, k-1)$
$K_1, class -: (x_0, y_0) = (k^2-k+1,1-k)$
$K_2, class +: (x_0, y_0) = (x_n, +y_n)$ for any $y_n \in K_1 \cup K_2$
$K_2, class -: (x_0, y_0) = (x_n, -y_n)$ " "
and in all cases above $(x_1, y_1)$ satisfy
$x_1 = Rx_0 + DSy_0$
$y_1 = Ry_0 + Sx_0$
Now if Andrej's conjecture is true, and we believe it is, then each operation "
add $y_n$" always adds a new $y_n$ to its list, and the two lists $K_1, K_2$ have
no common elements.
An implementation of
Enum_K with a bailout parameter finds
that with $k <10^6$ we have $|K_1| = 882, |K_2| = 163, |K| = 1045$, and
that every $k$ enumerated was unique. This agrees with Andrej's figure.
Best Answer
Just noticed this question. I agree with L.H.Gallardo that the problem is old (see e.g. Problem D5 in UPINT = Unsolved Problems in Number Theory by R.K.Guy), but not that it is hopeless: the usual heuristics suggest that the number of solutions with $\max(|x|,|y|,|z|) \leq H$ should be asymptotic to a multiple of $\log H$, so further solutions should eventually emerge (though it may indeed be hopeless to prove anything close to the $\log H$ heuristic).
See also my article
Rational points near curves and small nonzero $|x^3-y^2|$ via lattice reduction, Lecture Notes in Computer Science 1838 (proceedings of ANTS-4, 2000; W.Bosma, ed.), 33-63 = math.NT/0005139 on the arXiv.
Among other things it gives an algorithm for finding all solutions of $|x^3 + y^3 + z^3| \ll H$ with $\max(|x|,|y|,|z|) \leq H$ that should run (and in practice does run) in time $\widetilde{O}(H)$; since we expect the number of solutions to be asymptotically proportional to $H$, this means we find the solutions in little more time than it takes to write them down.
D.J.Bernstein has implemented the algorithm efficiently, and reports on the results of his and others' extensive computations at http://cr.yp.to/threecubes.html .
EDIT: for the specific problem $x^3+y^3+z^3=3$, Cassels showed that any solution must satisfy $x\equiv y\equiv z \bmod 9$ in this brief article:
A Note on the Diophantine Equation $x^3+y^3+z^3=3$, Math. of Computation 44 #169 (Jan.1985), 265-266.
This uses cubic reciprocity, and is stronger than what one can obtain from congruence conditions. See also Heath-Brown's paper "The Density of Zeros of Forms for which Weak Approximation Fails" (Math. of Computation 59 #200 (Oct.1992), 613-623), where he gives corresponding conditions for the homogeneous equation $x^3 + y^3 + z^3 = 3w^3$ and also $x^3 + y^3 + z^3 = 2w^3$, and reports that
In a letter to the author, Professor Colliot-Thélène has shown that the above congruence restrictions are exactly those implied by the Brauer-Manin obstruction. Moreover, for the general equation $x^3 + y^3 + z^3 = kw^3$, with a noncube integer $k$, there is always a nontrivial obstruction, eliminating two-thirds of the adèlic points.