Solution of Differential Equation by Numerical method in MATLAB

I am providing a matlab code, please help me why the displacement vs time plot gives the displacement increasing with time. though the excitation focre applied is sinusoidal, it was expected that it will be in oscillatory motion.

code:

% System parameters
    params.m = 0.91;               % Inner mass (kg)
    params.M = 1.2;                % Outer mass (kg)
    params.K = 250;            % Spring stiffness (N/m)
    params.c_m = 0.1;              % Increased mechanical damping coefficient (N·s/m)
    params.mu = 3.55e-7;           % Magnetic dipole moment (A·m^2)
    params.a = 0.0171;             % Coil radius (m)
    params.R = 0.121;              % Total resistance (Ohms)
    params.b0 = 0.005;             % Initial offset (m)
    params.F0 = 50;                % Reduced amplitude of external force (N)
    params.omega = 32;             % Angular frequency (rad/s)
    % Initial conditions [x, dx, y, dy]
    Z0 = [0, 0, 0.01, 0];          % Reduced initial displacement
    % Time span for simulation
    tspan = [0, 2];               % Simulation time range (10 seconds)
    % Solver options
    options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
    % Solve using ode15s
    [T, Z] = ode15s(@(t, Z) dynamic_equations(t, Z, params), tspan, Z0, options);
    % Extract results
    x = Z(:, 1);                   % Inner mass displacement
    dx = Z(:, 2);                  % Inner mass velocity
    y = Z(:, 3);                   % Outer mass displacement
    dy = Z(:, 4);                  % Outer mass velocity
    % Compute relative displacement, velocity, voltage, power, and forces
    z = params.b0 + (x - y);       % Total magnet-to-coil distance
    dz = dx - dy;                  % Relative velocity
    V = compute_voltage(z, dz, params);  % Voltage
    P = V.^2 / (params.R / 2);     % Power
    restoring_force = params.K * (x - y); % Spring restoring force
    external_force = params.F0 * sin(params.omega * T); % External harmonic force
    % Plot all results in a single figure
    plot_all_results_single_figure(T, x, dx, y, dy, z, V, P, restoring_force, external_force);
end
function dZdt = dynamic_equations(t, Z, params)
    x = Z(1); dx = Z(2); y = Z(3); dy = Z(4);
    % Total distance and velocity
    z = params.b0 + (x - y);
    dz = dx - dy;
    % Electromagnetic damping force
    C_e = (36 * pi^2 * params.mu^2 * params.a^2 * z^2) / (params.R * (params.a^2 + z^2)^(5/2));
    F_em = C_e * dz;
    % Accelerations
    ddx = (-params.K * (x - y) - params.c_m * dz - F_em) / params.m;
    ddy = (params.K * (x - y) + params.c_m * dz + F_em + params.F0 * sin(params.omega * t)) / params.M;
    % Return derivatives
    dZdt = [dx; ddx; dy; ddy];
end
function V = compute_voltage(z, dz, params)
    V = (6 * pi * params.mu * params.a^2 * z .* dz) ./ (params.a^2 + z.^2).^(5/2);
end
function plot_all_results_single_figure(T, x, dx, y, dy, z, V, P, restoring_force, external_force)
    % Create a single figure with six subplots
    figure;
    % Subplot 1: Displacement
    subplot(3, 2, 1);
    plot(T, x, 'r', 'LineWidth', 1.5); hold on;
    plot(T, y, 'b--', 'LineWidth', 1.5);
    xlabel('Time (s)');
    ylabel('Displacement (m)');
    legend('Inner Mass (x)', 'Outer Mass (y)');
    title('Time vs Displacement');
    grid on;
    % Subplot 2: Velocity
    subplot(3, 2, 2);
    plot(T, dx, 'r', 'LineWidth', 1.5); hold on;
    plot(T, dy, 'b--', 'LineWidth', 1.5);
    xlabel('Time (s)');
    ylabel('Velocity (m/s)');
    legend('Inner Mass (dx)', 'Outer Mass (dy)');
    title('Time vs Velocity');
    grid on;
    % Subplot 3: Relative Displacement (z)
    subplot(3, 2, 3);
    plot(T, z, 'g', 'LineWidth', 1.5);
    xlabel('Time (s)');
    ylabel('Relative Distance z (m)');
    title('Time vs Relative Distance');
    grid on;
    % Subplot 4: Voltage
    subplot(3, 2, 4);
    plot(T, V, 'k', 'LineWidth', 1.5);
    xlabel('Time (s)');
    ylabel('Voltage (V)');
    title('Time vs Voltage');
    grid on;
    % Subplot 5: Power
    subplot(3, 2, 5);
    plot(T, P, 'm', 'LineWidth', 1.5);
    xlabel('Time (s)');
    ylabel('Power (W)');
    title('Time vs Power');
    grid on;
    % Subplot 6: Forces
    subplot(3, 2, 6);
    plot(T, restoring_force, 'g', 'LineWidth', 1.5); hold on;
    plot(T, external_force, 'm', 'LineWidth', 1.5);
    xlabel('Time (s)');
    ylabel('Force (N)');
    legend('Restoring Force (K(x-y))', 'External Force (F_0 sin(\omega t))');
    title('Forces Comparison');
    grid on;
end

[Walter Roberson]

You have an extra end after the call to plot_all_results_single_figure

[Walter Roberson]

    ddy = (params.K * (x - y) + params.c_m * dz + F_em + params.F0 * sin(params.omega * t)) / params.M;

That line contains the only oscillary force, in the sin() term. But there are other forces being added to it, and those other forces have the potential to be larger than the sin() term.