mu=h^2
clear x y
x(1)=x0
y(1)=y0
plst='.r'
%  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;
  yh=y(k-1)-mu*sin(x(k-1))/2;
  for i=1:10
   % xo=xh
     xh =x(k-1)+yh/2;
     yh=y(k-1)-mu*sin(xh)/2;
   %  if abs(xo-xh)< 1e-10
  end      
  y(k)=y(k-1)-mu*sin(xh);
  x(k)=x(k-1)+yh;
 % 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
% y(k)=y(k)-(2*pi)*floor(y(k)/(2*pi));
% if(y(k)>pi)
%    y(k)=y(k)-2*pi;
%  end
end
plot(x,y,plst)
plot(-x,-y,plst)

