% phase plot of pendulum

% damping
B = .6
f = @(t,x) [x(2); -sin(x(1)) - B*x(2)];

y1 = linspace(-3,10,20);
y2 = linspace(-3,3,20);
[x,y] = meshgrid(y1,y2);

u = zeros(size(x));
v = zeros(size(x));

t=0;
for i = 1:numel(x)
    dx = f(t,[x(i); y(i)]);
    u(i) = dx(1);
    v(i) = dx(2);
end

quiver(x,y,u,v,'r'); figure(gcf)
xlabel('x_1')
ylabel('x_2')
axis tight equal;

hold on
for x = -2:2:8
  for y = -2:0.5:2
    [ts,ys] = ode45(f,[0,50],[x;y]);
    plot(ys(:,1),ys(:,2))
  end
end
hold off