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.1;0.];

options = odeset('RelTol',1e-6,'AbsTol',1e-6,'MaxStep',0.01);
[t,x]=ode45(@harmf,[0,60],x0,options);
%plot(t,x(:,1),'o')
%plot(t,[ekokku,epot,ekin]);
%legend('kokku','pot','kin');
%xx=1.*sin(x(:,1));yy=-1.*cos(x(:,1));
%xx=x(1)*cos(omega*t);

%plot(t,[x(:,1),xx(:)])
%egend('1','2');
%axis equal
%comet(xx,yy)
grid on;
plot(t,x(:,1))
hold off;
%a=load('out.dat');
%plot(a(:,1),a(:,2))
%hold off
%xlabel('time(s)');ylabel('energies');