


function dy = rhs(t,y)
global k;
m1 = 6;
m2 = 6;
dy = zeros(4,1);    % a column vector
dy(1) = (-k*y(3) + k*y(4) + Fext)/m1;
dy(2) =  (k*y(3) - k*y(4) - Fext)/m2;
dy(3) =  y(1);
dy(4) =  y(2);
end



function Fext = Fext()
Fext = 1;  %six lbs buoyant
end
