fdpend.m

function pos = fdpend(interval, ic, n)
h = (interval(2) - interval(1))/n;

y(1,:) = ic;
t(1) = interval(1);
direction = [0];
global pos
for k = 1:n
    t(k + 1) = t(k) + h;
    y(k + 1, :) = trapstep( t(k), y(k,:), h);
    
    if y(k+1,2)*y(k,2) < 0
        direction = [direction, y(k+1,1)];
    end

    
end

direction = direction(end:-1:1);
pos = [direction(1), direction(2)];
function y = trapstep(t,x,h)
% one step of the Trapezoid Method
z1 = ydot(t, x);
g = x + h*z1;
z2 = ydot(t + h, g);
y = x + h*(z1 + z2)/2;

function z = ydot(t,y)
d = 1;
A = 10;
g = 9.81; length = 1;
z(1) = y(2);
z(2) = -(g/length)*sin( y(1) ) - d*y(2) + A*sin(t);

samplingpend.m

data_vec = fdpend([0, 60], [0, 0], 12000);
for k = 1:24
    
    data_vec = [data_vec; fdpend([0, 60], [k*pi/12, 0], 12000)];
end


for k=1:25
    if data_vec(k, 1) > 2*pi
        while data_vec(k, 1) > 2*pi
            data_vec(k, 1) = data_vec(k, 1) - 2*pi;
        end
    elseif data_vec(k, 1) < 0
        while data_vec(k, 1) < 0
            data_vec(k, 1) = data_vec(k, 1) + 2*pi;
        end
    end

    if data_vec(k, 2) > 2*pi
        while data_vec(k, 2) > 2*pi
            data_vec(k, 2) = data_vec(k, 2) - 2*pi;
        end
    elseif data_vec(k, 2) < 0
        while data_vec(k, 2) < 0
            data_vec(k, 2) = data_vec(k, 2) + 2*pi;
        end
    end
end

data_vec