global omega2 beta;

k=25;m=1.;alpha=0.1;
omega2=k/m;
beta=2.*alpha/m;

x0=[0.;1.];

options = odeset('RelTol',1e-9,'AbsTol',1e-9,'MaxStep',0.001);
[t,x]=ode45(@harm,[0,3],x0,options);


subplot(1,2,1)
plot(t,x(:,1))
grid on;
subplot(1,2,2)
plot(x(:,1),x(:,2))
grid on;

function dx=harm( t,x )
global omega2 beta;
dx=zeros(2,1);
dx(1)=x(2);
dx(2)=-omega2*x(1)-2.*beta*x(2)/m;
end

