Find all ordered pairs $(x,y)$ of positive integers such that the expression $x^2+y^2+xy$ is a perfect square. My approach so far:
Let $x^2+y^2+xy=n^2,$ where $n\in\mathbb Z$.
$\implies (x+y) ^2-xy=n^2$
$\implies (x+y) ^2-n^2=xy$
$\implies (x+y+n) (x+y-n) =xy$
Only one case is possible:
When $x+y+n=xy$ and $x+y-n=1$. On adding, we get: $2(x+y) =xy+1$. Since LHS is an even number so is RHS. That's why, $xy$ must be odd$\implies x$ and $y$ both are odds. So I checked for few odd numbers and found out that for $x=3$ and $y=5,\;x^2+y^2+xy$ becomes $49$ which is a perfect square.. I didn't dare to check for more odds as there could be many...
Recently I found out that for $x=5$ and $y=16$(an even number), $x^2+y^2+xy=361=19^2.$ (Surprising!!) 
So now I can say, I am stuck very badly.. All of my observations miserably failed...
Please suggest something! 

Firstly, we can ssume Multiplying by $4$, we get $(2x+y)^2+3y^2$ is a perfect square. Thus we want to find $a,b$ such that $a^2+3b^2$ is a perfect square with $gcd(a,b)=1$ 
What follows is a well-known theory of finding out rational points on conic sections.
Now this amounts to finding rational points on the ellipse $E: X^2+3Y^2=1$. 
We know $(1,0)$ is a rational point. If $P$ is another such point, then the slope of the line joining $P$ and $(1,0)$ is rational $m$. Say $L$ intersects $E$ at $(u,v)$. Then we get $$ \frac{v}{u-1}=m
$$ and $$u^2+3v^2=1$$
This gives us $u^2+3m^2(u-1)^2=1$
Since $1$ ia anyway a root of the quadratic, the other root is $\frac{3m^2-1}{3m^2+1}$ and this gives $u=\frac{3m^2-1}{3m^2+1},v=\frac{-2m}{3m^2+1}$ 
Thus all rational points are parametrized by $\left  (\frac{3m^2-1}{3m^2+1},\frac{-2m}{3m^2+1} \right ) \ ; \ m\in \mathbb Q$ or $\left  (\frac{3m^2-1}{3m^2+1},\frac{2m}{3m^2+1} \right ) \ ; \ m\in \mathbb Q$
So we get if $a^2+3b^2=c^2$, then $\frac{a}{c}=\frac{3m^2-1}{3m^2+1}$ $\frac{b}{c}=\frac{2m}{3m^2+1}$ for some $m\in \mathbb Q$
writing everything in terms of integers we get $$a=\pm(3p^2-q^2),b=\pm 2pq, c= \pm (3p^2+q^2)$$ 
Thus we get $y=6kpq$ where $p,q,k>0$ and $2x+y=k|3p^2-q^2|$ 
If $p,q$ have opposite parity, this forces $k$ to be even and we get $k=2l$ 
So we finally get $y=4lpq$ and $x=l|3p^2-q^2|-2lpq$ 
If, on the other hand both $p,q$ are odd then we get $y=2kpq, x=\frac{k}{2}|3p^2-q^2|-kpq$ 
