mu=h^2
clear x y
x(1)=x0
y(1)=y0
% hold off
plot(x0,y0,plst)
axis([-pi pi -pi pi])
hold on
for k=2:np
  xh=x(k-1)+y(k-1)/2;
  y(k)=y(k-1)-mu*sin(xh);
  x(k)=xh+y(k)/2;
% riduzione al toro
  x(k)=x(k)-(2*pi)*floor(x(k)/(2*pi));
  if(x(k)>pi)
     x(k)=x(k)-2*pi;
  end
end
plot(x,y,plst)
% sfruttando la simmetria della mappa raddoppio il numero di punti
plot(-x,-y,plst)
