What numerical methods can solve $\sin(x) + \sin(y) = \sin(xy)$ Here is a nice graph representing the solution: wolframalpha.
I wish to draw such a graph myself but don't have any idea which methods exist and which of them are more appropriate for equations of this kind.
 A: Note that the plot from Wolfram Alpha isn't particularly good; you can see a lot of jagged artifacts all over the place, presumably due to insufficient initial sampling.  Here's a better one from Maple, made using

plots:-implicitplot(sin(x)+sin(y)-sin(x*y), x=-5..15, y=-10..10, grid=[100,100], gridrefine=2);


A: The standard numerical method for plotting an implicit function graph, say $f(x,y) = 0$, is to evaluate $f$ at points on a triangular grid, see which cells of the grid have a change of sign, and use linear interpolation to approximate the function in each cell.  Thus if you have three grid points $A$, $B$, $C$ with $f(A) = 1$, $f(B) = -1$ and $f(C) = 2$, you would draw a line segment from the midpoint of $AB$ to the point $1/3$ of the way from $B$ to $C$.
More sophisticated software might sample more points in a cell once a sign change has been detected, or use quadratic interpolation.
Another approach would be to numerically solve a system of differential equations such as $\frac{dx}{dt} = \frac{\partial f}{\partial y}$, $\frac{dy}{dt} = - \frac{\partial f}{\partial x}$, starting at an initial solution.  This would run into trouble at points where both partial derivatives are $0$ (which typically happens at self-intersections of the curve).  Also, it's hard to know a priori how many initial solutions to look for (one in each component of the curve).
A: The figure shows that there is an arc through the origin which near $(0,0)$ is well approximated by the hyperbola $x+y=x y$. For this arc (a tiny portion of the whole image!) you can recursively compute a Taylor expansion as follows: Write
$$x(t):=t +\sum_{k=0}^\infty a_{2k} t^{2k}\ ,\qquad y(t):=x(-t)$$
(note that I have chosen a particular parametrization) and plug this "Ansatz" into the equation $\sin\bigl(x(t)\bigr)+\sin\bigl(y(t)\bigr)=\sin\bigl(x(t)y(t)\bigr)$. Using the Taylor expansion of $\sin$ and comparing coefficients for all even powers of $t$ the $a_{2k}$ can now be determined one for one.  Doing this with Mathematica I have obtained the following first coefficients:
$$x(t)=t - {t^2\over2} - {t^4\over8} + {t^6\over12} - {319 t^8\over5760} - {541 t^{10}\over26880}+\ldots\quad .$$
Now make a parametric plot of $t\mapsto\bigl(x(t),x(-t)\bigr)$ for $|t|\leq 1$, say, and you will see the arc in question.
I have told this story because it gives an "analytical" approach to such problems and is radically different from the ideas brought forward in the other comments and answers.
