function ex33ch2% Define the physical constants:global m1 m2 L1 L2 gm1 = 937.5;m2 = 312.5;L1 = 16;L2 = 16;g = 32;% First solve the linear model:options = odeset('Mass',@massL,'MassSingular','no','Events',@events);[tL,yL,teL,yeL,ieL] = ode45(@odesL,[0 2*pi],[-0.5; -1; 1; 2],options);if ~isempty(ieL)    fprintf('Linear model:    lower pallet lost at t = %6.4f.\n',...             teL(end));endplot(tL,yL(:,1),'o',tL,-0.5*cos(2*tL) - 0.5*sin(2*tL))legend('Computed','Analytical')title('Linear Model: \theta_1(t)')figureplot(tL,yL(:,3),'o',tL,cos(2*tL) + sin(2*tL))legend('Computed','Analytical',2)title('Linear Model: \theta_2(t)')% Now solve the nonlinear model:options = odeset('Mass',@massN,'MassSingular','no','Events',@events);[tN,yN,teN,yeN,ieN] = ode45(@odesN,[0 2*pi],[-0.5; -1; 1; 2],options);if ~isempty(ieN)    fprintf('Nonlinear model: lower pallet lost at t = %6.4f.\n',...             teN(end));endfigureplot(tL,yL(:,1),'b',tN,yN(:,1),'r')legend('Linear model','Nonlinear model')title('Computed \theta_1(t)')figureplot(tL,yL(:,3),'b',tN,yN(:,3),'r')legend('Linear model','Nonlinear model',2)title('Computed \theta_2(t)')%==================================================================function dydt = odesL(t,y)% Linear model:global m1 m2 L1 L2 gdydt = zeros(4,1);dydt(1) = y(2);dydt(2) = -(m1 + m2)*g*y(1);dydt(3) = y(4);dydt(4) = -m2*g*y(3);function A = massL(t,y)global m1 m2 L1 L2 gA = zeros(4);A(1,1) = 1;A(2,2) = (m1 + m2)*L1;A(2,4) = m2*L2;A(3,3) = 1;A(4,2) = m2*L1;A(4,4) = m2*L2;function dydt = odesN(t,y)% Nonlinear model:global m1 m2 L1 L2 gdydt = zeros(4,1);dydt(1) = y(2);dydt(2) = m2*L2*y(4)*y(4)*sin(y(3) - y(1)) - (m1 + m2)*g*sin(y(1));dydt(3) = y(4);dydt(4) = -m2*L1*y(2)*y(2)*sin(y(3) - y(1)) - m2*g*sin(y(3));function A = massN(t,y)global m1 m2 L1 L2 gA = zeros(4);A(1,1) = 1;A(2,2) = (m1 + m2)*L1;A(2,4) = m2*L2*cos(y(3) - y(1));A(3,3) = 1;A(4,2) = m2*L1*cos(y(3) - y(1));A(4,4) = m2*L2;%==================================================================function [value,isterminal,direction] = events(t,y)global m1 m2 L1 L2 gvalue = y(3) - pi/3;isterminal = 1;direction  = 0;