Matlab code for restoring force plot

%Plot for restoring Force of a simple pendulum
% Written by S. K. Dwivedy on 30th May 2012
% th= theta
%L= length of pendulum

L=1;
g=9.8;

th=-pi:pi/100:pi;
f=(g/L)*(th-(1/factorial(3))*th.^3); %upto cubic order
f3=(g/L)*sin(th); % Actual
f1=10*th; % linear approximation
f5=(g/L)*(th-(1/factorial(3))*th.^3+(1/120)*th.^5); %uto quintic order
f7=f5-(g/L)*(1/factorial(7))*th.^7; %uto 7th order

plot(th,f,th,f1, 'r' ,th,f3, 'v' ,th,f5, 'g' ,th,f7, 'b' )
grid on

xlabel( '\theta' )
ylabel( 'Restoring force' )

Fig. 2.1.3: Different approximation of the restoring force.