f = 0.05;
t_f = 10;
tspan = [0 t_f];
z0 = [0 ; 0];

options = odeset('stats','on');

[t, z] = ode23s(@(t,z) sys_13(t, z, f), tspan, z0, options);

figure(1)
set(gcf,'units','normalized','outerposition',[0 0 1 1])

subplot(2,1,1);
plot(t,z(:,1),'b');
title(['Time history of $y(t),,(ddot{u}_{gx}=-(2pi f)^2sin(2pi f t)) ,, f=, $' num2str(f) ' $ $'],'interpreter','latex','fontsize',20)
ylabel('$y$','interpreter','latex','fontsize',15)
xlabel('$t$','interpreter','latex','fontsize',15)
xlim(tspan)

grid on
hold on

subplot(2,1,2);
plot(t,z(:,2),'b')
ylabel('$dot{y}$','interpreter','latex','fontsize',15)
xlabel('$t$','interpreter','latex','fontsize',15)
xlim(tspan)

grid on
hold on

figure(2)
plot(t,z(:,1).*z(:,2),'b')

grid on
grid minor
hold on

function dz = sys_13(t, z , f)
dz = zeros(2,1);

m = 1;
k0 = 4;
% tol = 1e-6;
% if abs(z(1)*z(2)) > tol && z(1)*z(2) > 0
%     k1 = 4;
% else
%     k1 = 0;
% end

if z(1)*z(2) > 0
    k1 = 4;
else
    k1 = 0;
end

finput = ( -(2*pi*f)^2 )*sin( 2*pi*f*t );

dz(1) = z(2);
dz(2) = - finput - ((k0 + k1)/m)*z(1);
end