clear x y
x(1)=x0
y(1)=y0
%  hold off
plot(x0,y0,plst)
hold on
for k=2:np
  y1=y(k-1);
  x1=x(k-1);
  k1x=-sin(x1);
  x2=x1+h*y1/2;
  y2=y1+h*k1x/2;
  k2x=-sin(x2);
  x3=x1+h*y2/2;
  y3=y1+h*k2x/2;
  k3x=-sin(x3);
  x4=x1+h*y3;
  y4=y1+h*k3x;
  k4x=-sin(x4);
  x(k)=x(k-1)+h*(y1+y4+2*y2+2*y3)/6;
  y(k)=y(k-1)+h*(k1x+k4x+2*k2x+2*k3x)/6;
 % 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)
     % sfruttando al simmetria del sistema, raddoppio il numero di punti
plot(-x,-y,plst)






