One can start showing the following:

The integers $n$ which are of the form $x^2+xy+y^2$, for two relatively prime integers $x,y$ are precisely those positive integers occurring as divisors of $m^2+m+1$, for some integer $m$.

In other words, the polynomial $f(x)=x^2+x+1$ has the property that the positive divisors of the integers it represents are precisely those integers that can be properly represented by its homogenization $F(x,y)=x^2+xy+y^2$ ("properly" here refers to the condition $(x,y)=1$).

The proof uses the fact that the imaginary quadratic order $\mathbf{Z}[x]/f(x)$ has class number one, as anticipated by Elkies. It goes as follows:

Let $n$ be a positive divisor of $m^2+m+1$, for some integer $m$.
Consider the quadratic form in $x,y$ given by:

$Q(x,y)=\frac{m^2+m+1}{n}x^2-(2m+1)xy+ny^2$;

it has integer coefficients, positive definite, and has discriminant equal to $-3$ (in particular it is primitive).

Since $h(-3)=1$, there is only one reduced, positive definite quadratic form of discriminant $-3$. This is $E(x,y)=x^2+xy+y^2$. Therefore $Q$ and $E$ are properly equivalent (that is there is a determinant-one change of variables taking one into the other), and since $Q$ certainly properly represents $n$, so does $E$.

The converse is similar and uses the fact that if a quadratic form $Q$ properly represents an integer $n$, then $Q$ is properly equivalent to a form of the type $nx^2+bxy+cy^2$ (this is lemma 2.3 of Cox's wonderful book "Primes of the form x^2+ny^2").

Once you understand the positive integers $n$ that are properly represented by $E$, then you can get them all, after scaling by squares.

The original problem is then reduced to understanding those integers $n$ for which $x^2+x+1$ has a zero in $\mathbf{Z}/(n)$. Using the Chinese Remainder Theorem, this can be reduced to the case where $n$ is a prime power $p^s$. Then for $p\neq 3$ Hensel's Lemma tells you that your equation has a solution mod $p^s$ if and only if it has a solution mod $p$. With quadratic reciprocity you can conclude that any prime divisor of $n$ has to be congruent to $1$ mod $3$. I am at the moment missing how you can solve the equation $x^2+x+1=0$ in $\mathbf{Z}/3^s$, but I think that a version of Hensel's Lemma applies. I'll think about it.

[EDIT: The equation $x^2+x+1=0$ has a solution in $\mathbf{Z}/(3^s)$, with $s\geq 1$, if and only if $s=1$. (Therefore $x^2+xy+y^2$ represents properly only $3$ and $1$ as powers of $3$.) One can see this by checking that there is no solution for $s=2$, and therefore for $s>2$. More generally, if $p>2$ then $x^{p-1}+x^{p-2}+\ldots+x+1=0$ has no solutions in $\mathbf{Z}/(p^s)$ with $s>1$. The $p$--adic valuation of an integer of the form $x^{p-1}+x^{p-2}+\ldots+x+1=(x^p-1)/(x-1)$ is either zero or one.]

(Suggested reading. Cox's book quoted above and Serre's paper: $\Delta=b^2-4ac$)

congruent to $-1$ mod 4. For $x^2+xy+y^2$, likewise, the valuation must be even for all primes congruent to $-1 \bmod 3$. This works because the relevant quadratic imaginary rings ${\bf Z}[i]$ and ${\bf Z}[(1+\sqrt{-3})/2]$ have unique factorization. – Noam D. Elkies Oct 17 '11 at 17:18