2DOF Harmonic Oscillator - Is this right? (How would I do this with ODE45?)

I'm not entirely sure what's going on with this result. I don't think resonant frequency increases then decreases with depth.

I've attached my equations. If you know of an easier method for getting the same thing accomplished, please let me know.

I used this source as a reference for writing my script:

https://www.uml.edu/docs/Two-Degree-Vibration-Theory_tcm18-191227.pdf

clc; clear; close all;

for DEPTH = 0:0.01:0.05 % Buried depths
radius = 0.045;  % Object radius (m)
SA = pi * radius^2; % SURFACE AREA (m^2)
NAT_FREQ = 220; Q = 20; % NATURAL FREQUENCY (Hz) | QUALITY FACTOR

RHO_SOIL = 1540; % (kg/m^3)
C_SOIL = 274; att_L = 0; att_T = 0;
C_SOIL_L = C_SOIL; C1 = C_SOIL_L / (1 + att_L * 1i); % (m/s)
C_SOIL_T = 120; C2 = C_SOIL_T / (1 + att_T * 1i); % (m/s)

mm = 12;
Km = mm * (2 * pi * NAT_FREQ)^2;
Rm = (2 * pi * NAT_FREQ * mm) / Q;

% Assigns values based on depths.
if DEPTH == 0
    ms = 0; Ks2 = 0; Ks1 = 0;
else
    ms = (1/3) * (RHO_SOIL * (DEPTH * (pi * radius^2)))/SA;
    Ks2 = real(RHO_SOIL * C1^2 * SA / DEPTH);
    Ks1 = real(RHO_SOIL * C2^2 * SA / DEPTH); % (Pa/m)
end

t = 1; P = 1; F = P * SA;
MASS = [ms 0; 0 mm];
SPRING = [(Ks1 + Ks2) -Ks2; -Ks2 (Km + Ks2)];

% THIS IS WHERE THE ERROR MOST LIKELY OCCURS.
for FREQ = 1:1000
    OMEGA = 2 * pi * FREQ; OMEGA = transpose(OMEGA);
    INPUT = F * exp(1i * OMEGA * t); FORCE = [INPUT; 0];
    if (DEPTH == 0)
        Rs2 = 0; Rs1 = 0;
    else
        Rs2 = ((-1 * imag(Ks1) / OMEGA))/SA;
        Rs1 = -imag(Ks2) / OMEGA;
    end
    DAMPING = [(Rs1 + Rs2) -Rs2; -Rs2 (Rm + Rs2)];
    Z = [(-OMEGA^2 * MASS(1,1) + 1i * OMEGA * DAMPING(1,1) + SPRING(1,1)) (-OMEGA^2 * MASS(1,2) + 1i * OMEGA * DAMPING(1,2) + SPRING(1,2));
        (-OMEGA^2 * MASS(1,2) + 1i * OMEGA * DAMPING(1,2) + SPRING(1,2)) (-OMEGA^2 * MASS(2,2) + 1i * OMEGA * DAMPING(2,2) + SPRING(2,2))];
    H = inv(Z);
    xs(FREQ) = H(1,1) * FORCE(1) + H(1,2) * FORCE(2);
    if DEPTH == 0
        xm(FREQ) = H(2,2) * FORCE(1);
    else
        xm(FREQ) = H(2,1) * FORCE(1) + H(2,2) * FORCE(2);
    end
    vm(FREQ) = (FORCE(1) / SA)/Z(2,2);
    vs(FREQ) = (FORCE(1) / SA)/Z(1,1);
end

xs = transpose(xs); vs = transpose(vs);
xm = transpose(xm); vm = transpose(vm);

if DEPTH == 0
    plot(1:1000, abs(vm))
else
    plot(1:1000, abs(vm))
end
xlabel('Frequency (Hz)'); ylabel('Velocity (m/s)');
legend('0m', '0.01m', '0.02m', '0.03m', '0.04m', '0.05m')
hold on
end

Answer 1

[darova]

What about ode45?

% x(1) == xs
% x(2) == dxs
% x(3) == xm
% x(4) == dxm
% fun = @(t,x) [ Vs; d2xs; Vm; d2xm ];
fun = @(t,x) [x(2); (P0-Rs1*x(2)-Rs2*(x(2)-x(4))-Ks1*x(1)-Ks2*(x(1)-x(3)) )/ms; ...
              x(4); ( 0- Rm*x(4)+Rs2*(x(2)-x(4)) -Km*x(3)+Ks2*(x(1)-x(3)) )/mm ];
[t,x] = ode45(fun,[0 5], [0 0 0 0]);
xs = x(:,1);
plot(t,xs)