How can we solve the diophantine equation? How can we find all the primitive solutions of the diophantine equation $x^2+3y^2=z^2$ ? 
Some solutions are $(\pm n , 0 , \pm n)$ for $n\in \mathbb{N}$. How can we find also the other ones?

We can follow the stereographic projection method used to find Pythagorean triples. We start by dividing through by $z^2$, giving us the equation
$$\left(\frac{x}{z}\right)^2 + 3\left(\frac{y}{z}\right)^2 = 1.$$
Or, in other words, we search for rational points on the ellipse
$$x^2 + 3y^2 = 1.$$
We'll take the point $(1, 0)$ and stereographically project onto the $y$-axis. Consider an arbitrary point $(0, y)$ on the $y$-axis, and consider the line
$$r = (1, 0) + t(-1, y)$$
between $(1, 0)$ and $(0, y)$. This line passes through the ellipse at $(1, 0)$ and a second point, which will be the point that stereographically projects onto $(0, y)$. The $t$ value for this point must satisfy,
$$(1 - t)^2 + 3(ty)^2 = 1 \iff t((3y^2 + 1)t - 2) = 0.$$
The $t = 0$ solution produces $(1, 0)$, so we discard it. The other solution is
$$t = \frac{2}{3y^2 + 1},$$
which yields the point on the ellipse,
$$\left(\frac{3y^2 - 1}{3y^2 + 1}, \frac{2y}{3y^2 + 1}\right).$$
Now, when we have a rational point on the ellipse, the line will be of rational slope, and hence will stereographically project to a point $(y, 0)$ where $y$ is rational. That is, the rational point on the ellipse must take the above form where $y \in \mathbb{Q}$. If we take $y = \frac{m}{n}$, where $n \neq 0$, and $m, n \in \mathbb{Z}$, then this becomes
$$\left(\frac{3\left(\frac{m}{n}\right)^2 - 1}{3\left(\frac{m}{n}\right)^2 + 1}, \frac{2\frac{m}{n}}{3\left(\frac{m}{n}\right)^2 + 1}\right) = \left(\frac{3m^2 - n^2}{3m^2 + n^2}, \frac{2mn}{3m^2 + n^2}\right).$$

To wrap this up, suppose we have an integer solution $(x, y, z)$ to the Diophantine equation $x^2 + 3y^2 = z^2$. If $z = 0$, then $x = y = 0$, which is one solution. Otherwise, $z \neq 0$ and $\left(\frac{x}{z}, \frac{y}{z}\right)$ must be a (rational) point on the ellipse. This point may be $(1, 0)$, in which case $(x, y, z) = k(1, 0, 1)$ (the points you've mentioned). The other solutions must take the form
$$\left(\frac{x}{z}, \frac{y}{z}\right) = \left(\frac{3m^2 - n^2}{3m^2 + n^2}, \frac{2mn}{3m^2 + n^2}\right).$$
Equating the fractions, we must have, for some integers $k, l$,
\begin{align*}
x &= k(3m^2 - n^2) \\
z &= k(3m^2 + n^2) \\
y &= l(2mn) \\
z &= l(3m^2 + n^2).
\end{align*}
As such, $k = l$, so our solution is of the form
$$(x, y, z) = k(3m^2 - n^2, 2mn, 3m^2 + n^2)$$
for integers $m, n, k$. Importantly, each of these are a solution, as
$$(k(3m^2 - n^2))^2 + 3(k(2mn))^2 = k^2(9m^4 + n^4 - 6m^2 n^2 + 12m^2 n^2) = (k(3m^2 + n^2))^2.$$
Note that the $(0, 0, 0)$ solution can be obtained from $k = 0$, but there is no way to obtain the $k(1, 0, 1)$ solutions. Hence, our general solution is
$$(x, y, z) = k(1, 0, 1) \text{ or } k(3m^2 - n^2, 2mn, 3m^2 + n^2) \text{ for } k, m, n \in \mathbb{Z}.$$
