


function dy = rhs(t,y)

t

rho = 1025;  %kg/m^3
g = 9.81;  %m/s^2
h = 100;  %depth (m)
press = 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 = 100;

%Piston-Screw
k_ps = 0.25*pi*d_s^2*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 = 0.25*pi*d_n^2*E_steel;
b_nh = 2*sqrt(k_nh*m_h);  % For Critical Damping

Fw_p = press*pi*d_p^2/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)) + Fw_h)/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(1)-v(2)) + Fw_p)/m_p;
dy(3) = (k_sn*(y(4)-y(3)) + b_sn*(v(4)-v(3)) ...
        +k_ps*(y(2)-y(3)) + b_ps*(v(2)-v(3)))/m_s;
dy(3) = (k_nh*(y(1)-y(4)) + b_nh*(v(3)-v(4)) ...
        +k_sn*(y(3)-y(4)) + b_sn*(v(3)-v(4)))/m_n;
dy(5:end) = v
end

