opspend.m

function y = opspend(interval, ic, n, P)
h = (interval(2) - interval(1))/n;



y(1,:) = ic;
t(1) = interval(1);



for k = 1:n
    t(k + 1) = t(k) + h;
    y(k + 1, :) = trapstep( t(k), y(k,:), h, P);
end

function y = trapstep(t,x,h, P)
% one step of the Trapezoid Method
z1 = ydot(t, x, P);
g = x + h*z1;
z2 = ydot(t + h, g, P);
y = x + h*(z1 + z2)/2;

function z = ydot(t,y, P)
d = 0.1;
A = P; %25
g = 9.81; length = 2.5;
z(1) = y(2);
z(2) = -( (g/length) + A * cos( 2*pi*t ) )*sin(y(1)) - d*y(2);

stable.m

format long
a = 12;
b = 13;
x = 0;
for steps = 1:100
    param = (a+b)/2;
    data = opspend([0, 60], [3.1, 0], 12000, param);

    length = size(data, 1);

    for k = 1:length
        if or(data(k, 1) > pi/2, data(k, 1) < -pi/2)
           x = 1;
        else
            x = 0;
        end
    end
    if x == 1
        b = (a + b)/2;
    else
        a = (a + b)/2;
    end
        
end
a
b