@Khang Ngo, I get the same graph as you, when I run your code.
Your system is equivalent to a mass-dashpot-spring system in which m=1, b=b, and k=w0^2. Since w0=10, k=100. The theoretical solution is
where yss(t)=steady state response to the driving force, and yh(t)= homogeneous solution =response to the initial conditions.
The homogeneous solution, yh(t), is the same in all four cases considered here, because the only thing that changes is the frequency of the inhomogeneous driving force. The homogeneous solution is
where 
and 
=10.000 and, from the initial conditions, C=5 and D=0. So we have
and =10.000 and, from the initial conditions, C=5 and D=0. So we have
.The steady state solution for driving force 
, when ω>0, is
, when ω>0, is 
where 
and 
. I obtained that result by standard calculus, as you can verify.
and . I obtained that result by standard calculus, as you can verify. In the special case where 
: B=0, 
, and therefore 
simplifies to 
.
: B=0, , and therefore simplifies to .The theoretical solution above and the Matlab symbolic solution are different, except in the special case where w=0. I also computed the result by numerical integration, using ode45(). The numerical integration result matches the theoretical solution and does not match the symbolic solution. I conclude the symbolic solution is not working right in this problem. I don't know why.
