clear
clf
global omega02 alpha;

l=1.;g=9.814; alpha=0.5;

omega02=g/l;
x0=[0.;10];

options = odeset('RelTol',1e-6,'AbsTol',1e-6,'MaxStep',0.01);
[t,x]=ode45(@harmf,[0,60],x0,options);

grid on;
plot(t,x(:,1)/(2*pi))
hold off;

function dx=harmf( t,x )
global omega02 alpha;
dx=zeros(2,1);
dx(1)=x(2);
dx(2)=-omega02*sin(x(1))-alpha*x(2);
end