Looks like the input signal (green) is \[-\sin(t)\]; therefore, the input angular frequency is \[\omega=1\text{rad/s}\].
Have you tried the Laplace transform method. The Laplace transform of the input \[-\sin(t)\] is \[-1/(s^2+1)\].
The Laplace transform of the output is the product of the input and system transforms.
\[
Y(s) = \frac{-1}{s^2+1} \, \frac{1}{s^2+s}=\frac{-1}{s(s+1)(s^2+1)}=\frac{A}{s}+\frac{B}{s+1}+\frac{Cs+D}{s^2+1}
\]
where A, B, C and D may be determined by partial fraction expansion. The time domain version of the output is of the form
\[
y(t) = A u(t) + B e^{-t} u(t)+ (\text{sinusoid \, at \, 1 rad/s}) \, u(t)
\]
so the output will definitely have a constant due to step magnitude A and and exponential decay B and a sinusoid. Does this make sense? Try for yourself and see if you can figure out the complete solution. But I agree, the solution that is simulated looks reasonable to me.
v_c