pend.m

function pend(interval, ic, n)
h = (interval(2) - interval(1))/n;
figure(1)
videostring = sprintf('ic = [%.4f,%.4f], h = %.3f.mp4',ic(1),ic(2),h)
writerObj = VideoWriter(videostring, 'MPEG-4');
open(writerObj);


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

set(gca, 'xlim', [-1.2 1.2], 'ylim', [-1.2 1.2], ...
    'XTick', [-1 0 1], 'YTick', [-1 0 1], ...
    'Drawmode', 'fast', 'Visible', 'on', 'NextPlot', 'add');
set(gcf, 'Renderer', 'zbuffer');
cla;
axis square % make aspect ratio 1-1

bob = line('color', 'r', 'Marker', '.', 'markersize', 40, ...
    'erase', 'xor', 'xdata', [], 'ydata', []);

rod = line('color', 'b', 'LineStyle', '-', 'LineWidth', 3, ...
    'erase', 'xor', 'xdata', [], 'ydata', []);

for k = 1:n
    t(k + 1) = t(k) + h;
    y(k + 1, :) = trapstep( t(k), y(k,:), h);
    xbob = sin( y(k + 1, 1) ); ybob = -cos( y(k + 1) );
    xrod = [0 xbob]; yrod = [0 ybob];
    set(rod, 'xdata', xrod, 'ydata', yrod)
    set(bob, 'xdata', xbob, 'ydata', ybob)
    if and(mod(k,10) == 0, k > n/2)
        drawnow; pause(h)
        frame = getframe;
        writeVideo(writerObj, frame);
    end
end
close(writerObj);
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);