function func = ode45spen ( d, k );
% VisSPen
%
% This m-file is part of the material of the lecture
%   Scientific Computing and Software Tools
% at the Math Department of Halle University, Germany.
% Summer term 2004
%
% Author :      Prof. Dr. M. Arnold, arnold@mathematik.uni-halle.de
% Version of :  May 11, 2004
%
% Visualization of Stiff Pendulum time integration by ode45  
%
% Parameter
%   d      (input)  : damping constant of spring
%   k      (input)  : stiffness constant of spring
%
% Example
%   k = 1.0e+6;   d = sqrt ( k );
%   ode45spen ( d, k );

% -> time integration
alpha0 = 5 * pi/180;   
opt = odeset ( 'RelTol', 1.0e-5, 'AbsTol', 1.0e-5, 'Refine', 1 );
sol = ode45 ( @stiffpen, [ 0 10 ], [ alpha0, 0, 0, 0 ]', opt, d, k );

hmean = mean ( diff ( sol.x ) );

% -> eigenvalues
yast = [0 0 9.81/k 0]';
jac  = zeros ( 4, 4 );
del  = 1e-10;
for j=1:4, 
    y = yast;
    y(j) = y(j) + del;
    jac(:,j) = stiffpen ( 0, y, d, k )/del;
end;

ev = eig ( jac );

% -> plot dimensions
points =  100;
xmin   = -4;
xmax   =  1.3;
ymin   = -4.5;
ymax   =  4.5;
dyplot = (ymax-ymin)/75;

% -> plot
markdom ( explrk(2, 4, points+1, 0),      1, xmin, xmax, ymin, ymax, dyplot, 0 );
markdom ( explrk(2, 1, points+1, 1+4*i),  1, xmin, xmax, ymin, ymax, dyplot, 0 );
markdom ( explrk(2, 1, points+1, 1-4*i),  1, xmin, xmax, ymin, ymax, dyplot, 0 );
plot ( [xmin xmax], [0 0], '-k' );                    
plot ( [0 0], [ymin ymax], '-k' );     
plot ( hmean*real(ev), hmean*imag(ev), '+b', 'MarkerSize', 8 );
hold off;

axis ( [ xmin xmax ymin ymax ] );                     
set (gca, 'DataAspectRatio', [ 1 1 1 ] );                    
title ( sprintf ( 'ode45 for stiff pendulum,  k = %6.2e ,  d = %6.2e', k, d ) );   
xlabel ( 'Re z' );
ylabel ( 'Im z' ); 
