function func = solspen,
% SolSPen
%
% 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
%
% Plot program for the solution trajectories of stiff pendulum
%
% Parameters:
%   none
%
% Example:
%   solspen;                 % print -dpng ../png/solspen.png

%
% -> Steuerparameter Ausgabe Bildschirm / PostScript
graphwidth  = 28.0;
graphheight = 15.0;
fac         =  0.95;
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 ] );

% -> control parameters graphical output
tl    =  5.0;
lmax  =  0.015;
lmin  = -0.005;
almax = 10.0 * pi / 180;
almin = -almax;

% -> physical parameters
k = [ 1.0e3, 1.0e6 ];
d = [ 1.0e1, 1.0e2 ];

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

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

% -> plot
for iplot=1:2,

    % -> time integration
    sol = ode15s ( @stiffpen, [ 0 te ], y0, opt, d(iplot), k(iplot) );

    subplot ( 2, 2, iplot );
    yout = sol.y(3,:);   
    xcut = max ( sol.x(yout>lmax) );   if isempty(xcut),  xcut = 0;  end;
    yout(sol.x<=xcut) = NaN;
    plot ( sol.x, yout, 'b' );
    xlabel ( 't' );
    ylabel ( '\Delta_l' );
    title ( sprintf ( '        k = %6.1e,  d = %6.1e', k(iplot), d(iplot)  ) );
    axis ( [ 0 tl lmin lmax ] );

    subplot ( 2, 2, iplot+2 );
    plot ( sol.x, sol.y(1,:), 'g' );
    xlabel ( 't' );
    ylabel ( '\alpha' );
    title ( sprintf ( 'Steifes Pendel    k = %6.1e,  d = %6.1e', k(iplot), d(iplot)  ) );
    axis ( [ 0 tl almin almax ] );

end;   

%
% -> 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;
