clear
global omega02 omegav beta2 f0m;

k=25.;m=1.;omegav=10; alpha=0.1;f0=0.0;

omega02=k/m;beta2=alpha/m;f0m=f0/m;

x0=[0.;15];

options = odeset('RelTol',1e-6,'AbsTol',1e-6,'MaxStep',0.01);
[t,x]=ode45(@harmf,[0,60],x,options);
grid on;
plot(t,x(:,1))
hold off;
