function [time,th] = mypend(th1,th2,L,I,b,Del,T)

time = [0:Del:T];
th(1) = th1;
th(2) = th2;

for k=3:length(time)
    th(k) = 2*th(k-1) - th(k-2) - ...
        Del*b/I*(th(k-1)-th(k-2)) - ...
        Del^2*(9.81*L/I)*sin(th(k-1));
end

figure(1)
plot(time,th)
hold on