clear
global Omega alpha

alpha=2;
Omega=10;
theta0=acos(alpha);

figure(1)
tspan=[0 5];
y0=[.1;0];
[t,y]=ode45(@beadf,tspan,y0);
plot(t,y(:,1),'k')

figure(2)
[theta,omega]=meshgrid(0:.2:pi,-10:1:10);
v1=omega;
v2=Omega^2*(cos(theta)-alpha).*sin(theta);
quiver(theta,omega,v1,v2)
hold on
plot(y(:,1),y(:,2))
axis([0 pi -10 10])
xlabel('theta')
ylabel('omega')
hold off

figure(3)
[theta,omega]=meshgrid(-pi:.2:pi,-20:1:20);
v1=omega;
v2=Omega^2*(cos(theta)-alpha).*sin(theta);
quiver(theta,omega,v1,v2,'k')
hold on

for k=.1:.1:3
    y0=[k;0];
    [t,y]=ode45(@beadf,tspan,y0);
    plot(y(:,1),y(:,2),'k')
end

axis([-pi pi -20 20])
xlabel('\theta')
ylabel('\omega')
hold off

figure(4)
clf
hold on
tspan=[0 2.5];

y0=[.1;0];
[t,y]=ode45(@beadf,tspan,y0);
plot(t,y(:,1),'k')

y0=[.2;0];
[t,y]=ode45(@beadf,tspan,y0);
plot(t,y(:,1),'k--')

y0=[.3;0];
[t,y]=ode45(@beadf,tspan,y0);
plot(t,y(:,1),'k:')

y0=[.4;0];
[t,y]=ode45(@beadf,tspan,y0);
plot(t,y(:,1),'k-.')


plot(t,acos(alpha),'k--')

xlabel('t')
ylabel('\theta')
hold off
