2DOF System - What is wrong with the implemented differential equations here?

For the following system, you should see frequency and velocity decrease as depth is increased. However, my script is showing the opposite trend. I changed my Ks2 to have a negative height as a quick fix, but I don't know if that makes my results mathematically sound as far as the differential equations are concerned. Do you see an issue here?

RHOs = 1600; % Density (kg/m^3)
RADIUS = 0.1525; % Radius (m)
HEIGHT = 8; % Depth (in)
NAT_FREQ = 188; % Natural Frequency (Hz)

% Mass Parameters
Mm = 22.34;
Km = 3.11 * 10^7; 
Rm = 2730;

Cs = 265; % Wave Speed (m/s)
ATTENUATION = 0.06; % Attenuation

% Calculations-------------------------------------------------------------

AREA = pi * RADIUS^2; % Surface Area (m^2)
HEIGHT = HEIGHT * 0.0254; % Depth (m)

CsP = Cs / (1 + ATTENUATION * 1i); % Longitudinal Wave Speed (m/s)
CsS = 0.5 * CsP; % Transverse Wave Speed

K = CsP^2 * RHOs; % Bulk Modulus (N/m^2) - Using Longitudinal Wave Speed
G = CsS^2 * RHOs; % Shear Modulus (N/m^2) - Using Shear Wave Speed
FREQ = 1:1000; FREQ = transpose(FREQ); % Frequency Data (1 - 1000 Hz)
LAMBDA = CsS ./ FREQ; % Wavelength from 1 - 1000 Hz (m)
COMP = K.^-1; % Compressibility (m^2/N)

% Textbook Formulas
Ms = RHOs * AREA * HEIGHT; % Mass (kg)
Ks1 = ((8 * pi) ./ LAMBDA) * G * RADIUS; % Shear Spring 1 (N/m^2)
Ks2 = (1 / COMP) * (AREA / -HEIGHT); % Spring (m^3/N) [changed height to negative sign to make correct trend]

Ks2 = real(Ks2); Rs2 = imag(Ks2); % Parameters
Ks1 = real(Ks1); Rs1 = imag(Ks1); % Parameters

PRESSURE = 1; FORCE = PRESSURE * AREA; % Converts Pressure (Pa) to Force (N)
NAT_OMEGA = 2 * pi * NAT_FREQ; % Converts Natural Frequency to Angular Frequency (rad/s)

vs = []; vm = []; % Initializes arrays for the for-loop.

for j = 1:length(FREQ)
        OMEGA = 2 * pi * j;
        
        if (HEIGHT == 0) % NO DEPTH
            X = -1i * OMEGA * FORCE / (-Km + 1i * OMEGA * Rm + OMEGA^2 * Mm);
            vs = [vs X];
            
        else % WITH DEPTH
            Rs1 = 0; Ks1 = 0;
            A11 = Ks1 + Ks2 - 1i * OMEGA * (Rs1 + Rs2) - OMEGA^2 * Ms;
            A12 = Ks2 + 1i * OMEGA * Rs2;
            A21 = A12;
            A22 = Km + Ks2 - 1i * OMEGA * (Rs2 + Rm) - OMEGA^2 * Mm;
            A = [A11 A12; A21 A22];
            X = A \ [-1i * OMEGA * FORCE; 0];
            vs = [vs X(1)]; % Velocity
            vm = [vm X(2)]; % Velocity
    end
end


plot(1:1000, abs(vs)); xlim([0 500]);
title('Without Shear Parameters');
xlabel('Frequency (Hz)'); ylabel('Velocity (m/s/1Pa)');
% legend('0in', '2in', '4in', '8in', '16in')
hold on

[darova]

Doesn't look like you are trying to solve diff equations. Do you?

[onamaewa]

I'm using a laplace transform and then applying cramer's rule to solve the equations. A & X do those steps.

[Bjorn Gustavsson]

Why bother with Laplace transform, why not solve it analytically? It is a a system of linear ODEs...