function yp = stiffpen ( t, y, d, k ),
% StiffPen
%
% 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 :  Apr 28, 2004
%
% Right hand side for stiff pendulum
%
% Parameters:
%   t      (input)  : time
%   y      (input)  : actual positions and velocities
%   d      (input)  : damping constant of spring
%   k      (input)  : stiffness constant of spring
%   yp     (output) : actual value of the right hand side
%
% Example:
%   alpha0 = 5 * pi/180;   k = 1.0e3;   d = sqrt(k);
%   y0 = [ alpha0, 0, 0, 0 ]';
%   te = 10;
%   opt = odeset ( 'RelTol', 1.0e-5, 'AbsTol', 1.0e-5 );
%   sol = ode45 ( @stiffpen, [ 0 te ], y0, opt, d, k );
%   plot ( sol.x, sol.y);

% -> allocate memory for output
yp = zeros ( size(y) );

% -> physical parameters
g    = 9.81;
len0 = 1.0;

% -> read vector of unknowns
alpha  = y(1);
alphap = y(2);
dlen   = y(3);
dlenp  = y(4);

% -> evaluate right hand side
yp(1) = alphap;
yp(2) = - g/(len0+dlen) * sin(alpha);
yp(3) = dlenp;
yp(4) = (len0+dlen)*alphap + g*cos(alpha) - d*dlenp - k*dlen;