function ydot = twomassfun(t,y,m1,m2,k,b)
%
%  Function to return ydot for two-mass system
%  with lower mass possibly bottoming out
%
g = 9.81;
x1 = y(1);
x2 = y(2);
x1dot = y(3);
x2dot = y(4);
bottomed = 0;
if x2 <= 0
    f2 =                    %% total upward force in spring + damper
    x2 = 0;
    if f2 < m2*g
        x2dot = 0;
        bottomed = 1;
    end
end;
x1ddot =                    %% from eqn of motion
if bottomed
    x2ddot = 0;
else
    x2ddot =                %% from eqn of motion
end;
ydot = [x1dot x2dot x1ddot x2ddot]';