function itime = plotspen ( kspan ),
% PlotSPen
%
% 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 3, 2004
%
% Driver and plot program for stiff pendulum
%
% Parameters:
%   kspan  (input)  : vector of stiffness constant of spring
%   itime  (output) : matrix of cpu times for integration
%
% Example:
%   kspan = logspace ( 3, 7, 11 );
%   plotspen ( kspan );                        %   print -dpng ../png/plotspen.png

%
% -> Steuerparameter Ausgabe Bildschirm / PostScript
graphwidth  = 14.0;
graphheight = 15.0;
fac         =  0.9;
shift       =  0.5;
papershift  =  0.0;

ActFontSize   = 14;
ActMarkerSize =  6;
RatioFlag     =  0;

border      = 1.0;

%
% -> Bestimme Groesse des Bildschirms und Skalierungsfaktor fuer die Ausgabe
root = 0;
set ( root, 'Units', 'centimeters' );
screensize = get ( root, 'ScreenSize' );
left   = screensize(1);
bottom = screensize(2);
width  = screensize(3);
height = screensize(4);
set ( root, 'Units', 'pixels' );
if RatioFlag==0,
  screenscale = 1.0;
elseif RatioFlag==1,
  screenscale = (width-2*border) / graphwidth;
elseif RatioFlag==-1,
  screenscale = (height-2*border) / graphheight;
end;

%
% -> Definiere den Bildschirmausschnitt
if length ( get ( root, 'Children' ) ) == 0,
  h = figure;
else
  h = gcf;
  clf;
end;
stdunits = get ( h, 'Units' );
set ( h, 'Units',    'centimeters' );
set ( h, 'Position', ...
  [ left+width-screenscale*graphwidth-border height-2*border-screenscale*graphheight ...
    screenscale*graphwidth                   screenscale*graphheight ] );
set ( h, 'Units', stdunits );

%
% -> Skalierung der Ausgabe in PostScript
%
set ( gcf, 'PaperUnits',    'centimeters' );
set ( gcf, 'PaperPosition', [ papershift papershift graphwidth graphheight ] );

% -> get data
nk = length ( kspan );
itime = zeros ( nk, 2 );
tn    = cell ( nk, 2 );
h     = cell ( nk, 2 );

for ik=1:nk,

    % -> physical parameters
    k = kspan(ik);
    d = sqrt(k);

    % -> initial values
    alpha0 = 5 * pi/180;   
    y0 = [ alpha0, 0, 0, 0 ]';
    te = 10;

    % -> integrator settings
    opt = odeset ( 'RelTol', 1.0e-5, 'AbsTol', 1.0e-5, 'Refine', 1 );

    % -> time integration
    tic;
    sol = ode45 ( @stiffpen, [ 0 te ], y0, opt, d, k );
    tn(ik,1)    = {  sol.x(1:end-1) };
    h(ik,1)     = {  diff(sol.x) };
    itime(ik,1) = toc;
    
    tic;
    sol = ode15s ( @stiffpen, [ 0 te ], y0, opt, d, k );
    tn(ik,2)    = {  sol.x(1:end-1) };
    h(ik,2)     = { diff(sol.x) };
    itime(ik,2) = toc;

end;

subplot ( 2, 1, 1 );
semilogx ( kspan, itime(:,1), 'b', kspan, itime(:,2), 'g' );
legend ( 'ode45', 'ode15s', 2 );
xlabel ( 'k' );
ylabel ( 'cpu-time [s]' );
title ( 'Steifes Pendel' );
set ( gca, 'YLim', [ 0 3 ] );

subplot ( 2, 1, 2 );
semilogy ( tn{1,2}, h{1,2}, 'g', tn{1,1}, h{1,1}, 'b' );
legend ( 'ode15s', 'ode45', 4 );
hold on;
for ik=3:2:nk,
    semilogy ( tn{ik,2}, h{ik,2}, 'g', tn{ik,1}, h{ik,1}, 'b' );
end;
hold off;
xlabel ( 'Time  t' );
ylabel ( 'Stepsize  h' );
title ( 'Steifes Pendel' );
axis ( [ 0 te 1e-6 1 ] ); 
    
%
% -> Korrektur der Groesse und Lage der Koordinatensysteme, Aenderung der Schriftgroesse
allax = get ( gcf, 'Children' )';
for ax = allax,
  rect    = get ( ax, 'Position' );
  rect(2) = rect(2) + shift * (1-fac) * rect(4);
  rect(4) = fac * rect(4);
  set ( ax, 'Position', rect );
  set ( ax, 'FontSize', [ActFontSize] );
  set ( get ( ax, 'title' ), 'FontSize', [ActFontSize] );
  set ( get ( ax, 'xlabel' ), 'FontSize', [ActFontSize] );
  set ( get ( ax, 'ylabel' ), 'FontSize', [ActFontSize] );
  allline = get ( ax, 'Children' )';
  for al=allline,
    if get ( al, 'Type' ) == 'line',
      set ( al, 'MarkerSize', [ ActMarkerSize ] );
    end;
  end;
end;
