function mech_osc_demo(D,Method,pflag)
%
% Demonstrates MATLAB's ode45, ode113, and ode15s on the 
% mechanical oscillator IVP x'' = - Dx' - x, x(0) = x(1) = 1, 
% written as a first-order system 
%
%      y' = Ay, A = [0 1; -1 -D]
%      y_1(0) = y_2(0) = 1 
%
% Inputs: 
%  Required:
%    D      = drag coefficient. 
%  Optional: 
%    Method = 'ode45', 'ode113', or 'ode15s' (string). 
%    pflag  = plotting flag: 1 => solution vs time, 2 => phase plane. 
%
% The initial data are x(0) = y_1(0) = 1, x'(0) = y_2(0) = 1. 
%
% Interesting exercise: Compare the solution-vs-time and phase-plane 
% plots for D = 0, .1, 1, 2. Compare the performance of ode45, 
% ode113, and ode15s for D = 1, 10, 100, 1000.
%
% Required M-file: mech_osc_fun.m 
%

global D_GL; 

D_GL = D; 

% Set the defaults as needed. 

if nargin<3, pflag = 1; end 
if nargin < 2, Method = 'ode45'; end 

% Do the remaining initialization. 

y_0=[1;1]; 
tspan = [0 32];
options = odeset('Stats','on'); 

% Call the ODE solver. 

if strcmp(Method,'ode45')
  [T,Y] = ode45(@mech_osc_fun,tspan,y_0,options);
elseif strcmp(Method,'ode113')
  [T,Y] = ode113(@mech_osc_fun,tspan,y_0,options);
elseif strcmp(Method,'ode15s')
  [T,Y] = ode15s(@mech_osc_fun,tspan,y_0,options);
else
  error('Method must be ode45, ode113, or ode15s.');
end

% Plot the results. 

if pflag == 1, 
  plot(T,Y(:,1),T,Y(:,2)); 
  legend('x(t)','dx/dt(t)'); xlabel('t'); 
elseif pflag == 2
  plot(Y(:,1),Y(:,2));
  xlabel('x(t)'); ylabel('dx/dt(t)'); axis square;
else
  error('pflag must be 1 or 2.');
end 
title(['Damped Oscillator IVP,  D =',num2str(D),',  Method = ', Method]); 
figure(1); drawnow; 
