function bisection();
%% Find the initial angle theta0 such that the period is 10 seconds
mytol = 1e-6; % the tolerance
figure; grid; hold on;

%% initial bracket 
a = 0.99*pi; b = 2*pi/3;
[pa, Ta, Ya] = mypendulum(a); [pb, Tb, Yb] = mypendulum(b);
fa = pa-10; fb = pb - 10;
ha = plot(Ta(:,1),Ya(:,1),'r');
hb = plot(Tb(:,1),Yb(:,1),'k');
drawnow;
if sign(fa) == sign(fb)
    fprintf('Provided a and b do not bracket the solution\n');
    return;
end
%% while cycle
k=1;
while abs(b-a) > mytol
    c = 0.5*(b+a);
    fprintf('iter #%d,\t',k);
    if c==a | c==b
        fprintf('Solution is not found, mytol is too small\n');
        break;
    end
    [pc, Tc, Yc] = mypendulum(c);
    fc = pc - 10;
    fprintf(' period = %d\n', pc);
    if fc == 0
        fprintf('Solution is found: theta0 = %d, period = %d\n',c, fc+10); 
        break; % solution is found
    end
    if sign(fc) ~= sign(fb)
        a=c;
        fa=fc;
        plot(Tc(:,1),Yc(:,1),'r');
        drawnow;
        pause(1);
    else
        b=c;
        fb=fc;
        plot(Tc(:,1),Yc(:,1),'k');
        drawnow;
        pause(1);
    end
    k=k+1;
end
%% print the solution
if abs(b-a) <= mytol
    fprintf('Solution is found: theta0 = %d*pi, period = %d\n',c/pi, fc+10); 
end
%% Plot the resulting solution
options=odeset('AbsTol',1e-12,'RelTol',1e-12,'Events',@events);
f=@(t,x)[x(2),-0.98*sin(x(1))]';
[T,Y]=ode45(f,[0,50],[c,0],options);
plot(T(:),Y(:,1),'LineWidth',2);
% figure; grid; hold on;
% plot(Y(:,1),Y(:,2),'LineWidth',2);%

end
%% function period(theta0)
function [myperiod, T, Y]=mypendulum(theta0);
options = odeset('Events',@events,'AbsTol',1e-12,'RelTol',1e-12);
f = @(t,x)[x(2),-0.98*sin(x(1))]';
[T,Y] = ode45(f,[0,50],[theta0,0],options);
n = length(T);
myperiod = T(n);
end
%% stopping criterion: stop after one period
function  [value,isterminal,direction] = events(t,y)
value = y(2);
isterminal = 1; % stop integration if value = 0, i.e., if y(2) = 0
direction = -1; % stop integration only if the sign of y(2) switches from + to - 
end

