function dx=harmf( x )
global omega02 omegav beta2 f0m;
dx=zeros(2,1);
dx(1)=x(2);
dx(2)=-omega02*sin(x(1))-beta2*x(2)+f0m*cos(omegav*t);
end

