function dy = rhs(t,y)

rho = 1025;  %kg/m^3
g = 9.81;  %m/s^2
h = 100;  %depth (m)
p = rho*g*h;

E_steel = 195000000;  %Modulus of steel

m_s = .2;  %Screw Mass (kg)
m_h = 3;   %Housing Mass(kg)
m_p = .5;  %Piston Mass (kg)
m_n = .2;  %Nut Mass (kg)

d_s = .25*.0254;  %Screw Diameter  (m)
d_n = 1.0*.0254; %Nut Diameter (m)
d_p = 4*0.0254;  %Piston Diameter (m)



%Housing-Piston
k_hp = 0;
b_hp = .1;

%Piston-Screw
k_ps = d_s*E_steel;
b_ps = 2*sqrt(k_ps*m_s)  % For Critical Damping

%Screw-Nut
k_sn = 0.0;
b_sn = 1.0;

%Nut-Housing
k_nh = d_n*E_steel;
b_nh = 2*sqrt(k_nh*m_h);  % For Critical Damping

Fw_p = press*pi*d_p^4/4;
Fw_h = -Fw_p;

v = y(length(y)/2+1:end);

dy = zeros(length(y),1);    % a column vector
dy(1) = (k_hp*(y(2)-y(1)) + b_hp*(v(2)-v(1)) ...
        +k_nh*(y(4)-y(1)) + b_nh*(v(4)-v(1)))/m_h ;
dy(2) = (k_ps*(y(3)-y(2)) + b_ps*(v(3)-v(2)) ...
        +k_hp*(y(1)-y(2)) + b_hp*(v(2)-v(1)))/m_h ;
dy(3) =  y(1);
dy(4) =  y(2);
dy(5:end) = v;
end



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