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), direction(3), direction(4)];
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 = 12;
g = 9.81; length = 1;
z(1) = y(2);
z(2) = -(g/length)*sin( y(1) ) - d*y(2) + A*sin(t);

fdp12.m

marker = [];
storage = [];
for k = 1:24
    marker = [marker; k];
    
    data_vec_a = fdpend([0, 240], [k*pi/12, 0], 48000);

    for k=1:4
        if data_vec_a(1, k) > 2*pi
            while data_vec_a(1, k) > 2*pi
                data_vec_a(1, k) = data_vec_a(1, k) - 2*pi;
            end
        elseif data_vec_a(1, k) < 0
            while data_vec_a(1, k) < 0
                data_vec_a(1, k) = data_vec_a(1, k) + 2*pi;
            end
        end
    end
    storage = [storage; data_vec_a ];
end
res = [marker storage]

fdpab.m

marker = [];
storage = [];
for n = 1:24
    
    marker = [marker; n; n];
    
    data_vec_a = fdpend([0, 120], [n*pi/12 - 0.05, 0], 24000);
    data_vec_b = fdpend([0, 120], [n*pi/12 + 0.05, 0], 24000);

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

    for k=1:4
        if data_vec_b(1, k) > 2*pi
            while data_vec_b(1, k) > 2*pi
                data_vec_b(1, k) = data_vec_b(1, k) - 2*pi;
            end
        elseif data_vec_b(1, k) < 0
            while data_vec_b(1, k) < 0
                data_vec_b(1, k) = data_vec_b(1, k) + 2*pi;
            end
        end
    end
    storage = [storage; data_vec_a; data_vec_b];
end
results = [marker storage]