function duration
clear;
clc

theta=linspace(pi/360,350*pi/360,300);
T=[];
options=odeset('Events',@events);%开启事件判断功能
%解不同的初始角度下的周期值
for i=1:length(theta)
    [t,u]=ode45(@f,[0,50],[theta(i),0],options);
    T=[T,2*t(end)];
end
hold on
plot(theta,T,'r*-')
title('周期与摆角的关系');
xlabel('摆角'); ylabel('周期');

function ydot=f(t,y)
q=0; 
b=0; 
w0=2/3; %周期性强迫外力的频率
ydot=[y(2);-q*y(2)-sin(y(1))+b*cos(w0*t)];

function [value,isterminal,direction]=events(t,y)
value=y(2);   %监控的值为角速度
isterminal=1; %isterminal=1表示中断运算
direction=1;%direction=1表示由负变到零