function ydot = mbksys(t,y,m,b,k,Pfunc,pargs)
%MBKSYS -- solving single DOF 2nd order sys
%
%  Inputs coming from ode45
%       t               times
%       y               state vector (x, xdot)
%       m,b,k           mass, damping, stiffness
%       Pfunc(t,pargs)  function to get load P(t) (optional)
%       pargs           misc. inputs needed by Pfunc
%
%  Output
%       ydot      derivative of state vector
%
ydot = zeros(2,1);
if nargin > 5
    P = Pfunc(t,pargs);
else
    P = 0;
end
ydot(1) = y(2);
ydot(2) = (P - b*y(2) - k*y(1))/m;