Q:

How to solve an equation with two variables in MATLAB

For example, I want to solve equation |x|+|y|=1 numerically, I would to get a rectangle. If I plot it out, it should be something like 

It's a simple example. What I really want to solve is a equation like sin(x)sin(2y)+sin(y)sin(2x)=0.

A:

It really pays to do a little math before jumping to fancy graphs or numerical solvers.   
Suppose you have 
sin(x)·sin(2y) + sin(y)·sin(2x) = 0

From any standard list of trig identities you find that this can be re-written as
sin(x)·sin(y) · (cos(x) + cos(y)) = 0

So your equation holds when 
  sin(x)·sin(y) = 0   or 
cos(x) + cos(y) = 0

in other words,
x = 0  ± k·π,  or
y = 0  ± k·π,  or
x = -y ± 2k·π

with k ∈ ℤ0. 
Graphically, this would be a set of vertical lines (x = 0 ± k·π), a set of horizontal lines (y = 0 ± k·π) and a set of lines at ±45° (y = -x ± 2k·π).
Not a line of code needed. 

A:

With your example equation, Rody's answer provides an analytic solution.
In general, however, given an equation f(x,y)=0 it may not be possible to find an analytic solution. In that case I suggest:

Define a target x, y area and generate random samples within it.
Compute z = abs(f(x,y)) at the sample points.
Sort the values of z and define the approximate solution set as a given proportion of the lowest values.
f = @(x,y) sin(x).*sin(2*y)+sin(y).*sin(2*x); %// your example function
N = 1e7; %// number of samples
proportion = 1e-2; %// choose to achieve thin lines in solution set
xmin = -10; xmax = 10; %// x limits of target area
ymin = -10; ymax = 10; %// y limits of target area

x = xmin+(xmax-xmin)*rand(1,N);
y = ymin+(ymax-ymin)*rand(1,N);
z = abs(f(x,y));
[zsort isort] = sort(z);
n = round(N*proportion);
plot(x(isort(1:n)),y(isort(1:n)),'b.','markersize',.1)
axis square

