function my_ode_demo

%  Demonstrate the use of nested functions to pass information to MATLAB's
%  ODE solvers.
%    Solve the damped mass-spring model equation
%       m*x'' + gamma*x' + k*x = r(t)
%       x(0) = x0
%       x'(0) = v0
%  where external forcing function r(t) = A0*sin(omega_0*t).
%  System form is 
%               [  y2                             ]
%       yvec' = [                                 ]
%               [  -k/m*y1 - gamma/m*y2 - r(t)/m ]
%
%                  [  x0 ]
%       yvec(0) = [      ]
%                  [  v0 ]

%%  Set up parameters m=mass, gamma=damping_coef, k=spring_const.
  mass = 1;  
  damping_coef = 1;
  spring_const = 1;
  
%%  Set amplitude and frequency of forcing function 
%%    r(t) = A0*sin(omega_0*t).
  A0 = 1;
  omega_0 = 1;
  
%%  Set up reduced parameters k/m, gamma/m, and A0/m.
  kdm = spring_const/mass;
  gdm = damping_coef/mass;
  A0dm = A0/mass;
  
%%  Set initial conditions.
  x0 = 1;
  v0 = 0;
  yvec_init = [x0; v0];
  
%%  Set time interval.
  b = 40;
  time_interval = [0; b];
  
%%  Compute (constant) Jacobian.
  J = [0 1; -kdm -gdm];
  
%%  Set ODE solver options.
  options = odeset('Jacobian',J);
  
%%  Apply MATLAB ODE solver.
  [tvec,yvec] = ode15s( @rhs, time_interval, yvec_init, options );
  
%%  Display results
figure(1)
  plot(tvec,yvec(:,1),'-', tvec,yvec(:,2),'--')
  xlabel('time')
  ylabel('ODE solution components vs time')
  legend('displacement','velocity')
  
figure(2)
  plot(yvec(:,1),yvec(:,2))
  axis('equal')  %  insure proper aspect ratio
  xlabel('displacement')
  ylabel('velocity')
  title('phase portrait')
  
figure(3)
  plot(tvec(1:end-1),diff(tvec))
  xlabel('time')
  title('ODE solver step size vs time')
  
%% The following nested function evaluates the right-hand-side of the ODE
%% system. Parameter values are set "outside" in the main function.

    function fvec = rhs(t,yvec)
        
        y1 = yvec(1);
        y2 = yvec(2);
        fvec = zeros(2,1);  %  insure that fvec is a row and not a column vector
        fvec(1) = y2;
        fvec(2) = -kdm*y1 - gdm*y2 + forcing_fn(t);
        
    end

%%  This nested function evaluates the forcing function for the damped
%%  spring model.

    function r = forcing_fn(t)
        r = A0dm*sin(omega_0*t);
    end

end  %  end of main function