hold off;
c=0.3;
[x1,x2]=meshgrid(-2*pi:0.1:2*pi,-2:0.1:2);
px1=x2;
px2=sin(x1)-c*x2;
quiver(x1,x2,px1,px2);
f=@(t,X)[X(2); sin(X(1))-c*X(2)];
hold on;
for y0=[0 0.5 1 1.5 2 2.5]
    [ts,ys]=ode45(f,[0,50],[-pi;y0]);
    plot(ys(:,1),ys(:,2));
    plot(ys(1,1),ys(1,2),'bo');
    plot(ys(end,1),ys(end,2),'ks');
    [ts,ys]=ode45(f,[0,50],[pi;y0]);
    plot(ys(:,1),ys(:,2));
    plot(ys(1,1),ys(1,2),'bo');
    plot(ys(end,1),ys(end,2),'ks');
end
axis([-2*pi,2*pi,-2,2]);
