function unstable_soln_demo(theta_prime,Method,tspan)
%
% Demonstrates the instability of the solution of the pendulum 
% ODE 
%
%     theta'' + sin(theta) = 0, theta(0) = 0, theta'(0) = 2.
% 
% The integration uses ode113 with tight tolerances to show the 
% behavior of solutions in a narrow range of the unstable 
% solution above. 
%
% Inputs (optional): 
%  theta_prime = initial d(theta)/dt; may be a scalar or a vector. 
%  Method      = 'ode45' or 'ode113'. 
%  tspan       = integration range and output points. 
%
% Interesting exercises: First, call the routine with no arguments 
% and observe the solutions generated with the default values of 
% theta_prime. Then call the routine with theta_prime = 2 and 
% compare the performance of ode45 with that of ode113. Vary tspan 
% to see how far (in time) you can track the unstable solution with 
% ode113 and ode45. You can also change RelToland AbsTol in the 
% ``options = odeset'' line below to see how varying the tolerances 
% affects how far the unstable solution can be accurately tracked. 
% (They are currently set very small.) 
% 
% Required M-file: pendulum_fun.m 
%

% Set the defaults as needed. 

if nargin < 3, tspan = [0:.1:6*pi]; end 
if nargin < 2, Method = 'ode113'; end 
if nargin < 1, theta_prime = [1.98:.005:2.02]'; end 

% Complete the initialization. 

figure(1);clf;

options=odeset('RelTol',1.e-13,'AbsTol',1.e-14,'Stats','on');

% Loop over the theta_prime values. 

for i = 1:length(theta_prime);

% Call the ODE solver. 

  if strcmp(Method,'ode45'), 
    [T,Y]=ode45(@pendulum_fun,tspan,[0;theta_prime(i,1)],options);
  elseif strcmp(Method,'ode113'), 
    [T,Y]=ode113(@pendulum_fun,tspan,[0;theta_prime(i,1)],options);
  else
    error('Method must be ode45 or ode113.');
  end

% Plot the results. 

  if theta_prime(i,1)==2, 
    figure(1); plot(T,Y(:,1),'r-'); hold on; 
  else
    figure(1); plot(T,Y(:,1)); hold on; 
  end
end

title('Pendulum Problem Solutions'); drawnow; 
hold off; 

