t=0:0.001:200;
N=length(t);
dt=(t(end)-t(1))/N;
g=9.8; L=9.8; m=9.8;

A=0; B=0; gamma=0;
starting_angles=[pi/4, pi/2, .95*pi];

for k=1:length(starting_angles)
    theta(1)=starting_angles(k);
    omega(1)=0;
    KE(1)=(1/2)*m*L^2*omega(1)^2;
    PE(1)=m*g*((1-cos(theta(1)))*L);
    TE(1)=KE(1)+PE(1);

    for j=1:N-1
        theta_m=theta(j)+(omega(j)*dt);
        omega_m=omega(j) - (( (g/L * sin(theta(j))) + (gamma*omega(j)) - (A*cos(B*t(j))) )*dt);
        theta(j+1)=theta(j) + (omega_m*dt);
        omega(j+1)=omega(j) - (( (g/L * sin(theta(j))) + (gamma*omega(j)) - (A*cos(B*t(j))) )*dt);
        KE(j+1)=(1/2)*m*L^2*omega(j+1)^2;
        PE(j+1)=m*g*((1-cos(theta(j+1)))*L);
        TE(j+1)=KE(j+1)+PE(j+1);
    end

    figure(1)
    plot(theta,omega)
    hold on

    figure(2)
    plot(t,theta)
    hold on

    figure(3)
    plot(t,omega)
    hold on
end