% Consider the second order differential equation known as the
% pendulum equation
%             y'' + omegasq*sin(y) = 0 approximated as: y'' + omegasq*(y-(mu/6)*y^3)= 0
% where the prime ' denotes the time-derivative operator.  We can rewrite
% this as a system of first order differential equations:
%             y1' = y2
%             y2' = -omegasq*(y1-(mu/6)*y1^3)
%
% To simulate the above system, we create an m-file pendSim.m that calls
% a function M-file pendulum.m which returns the state derivative vector. It invokes ode23.m
% vdpsim.m plots teme responses of y and ydot in Figure 1 and the phase plane plot in Figure 2  
% It also queries new inputs for the parameter mu and the initial condition y(0)

clear all;

status='y';
mu=1.0; omegasq=16.;            % default values of parameters
tspan=[0.0 20.0];               % initial time t0 = 0; final time tfinal = 20;
y0 = [1.0 0.0]';                % default initial conditions y(0)=0  ydot(0)=0.1;

disp('    ')
disp('Pendulum Equation  d2y/dt2 + omegasq*(y-(mu/6)*y^3) = 0')
disp('    ')

while (status=='y')
   mu=input('enter the value of mu (default mu=1.0) ');
   y0(1)=input('enter the value of y0(1) (default y0(1)=1.0) ');
   options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4]);  %options=[];
   [t,y] = ode23('pendulum',tspan,y0,options,mu,omegasq);

   figure(1)
   plot(t,y) 
   title('Pendulum equation time history')
   xlabel('time t')
   ylabel('position y and velocity ydot')
   grid on
   
   figure(2)
   plot(y(:,1),y(:,2)) 
   title('Pendulum equation - phase plane plot')
   xlabel('position y')
   ylabel('velocity ydot')
   grid on
   
   status=input('enter y to continue or anything else to terminate  ','s')
end