function harm;
global omega;
omega=5.;
x0=[1.;0.];
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',1.);
[t,x]=ode45(@rpp,[0,200], x0,options);

plot(t,x(:,1))



function f=rpp(t,x)
global omega;
f=[x(2);-x(1)*omega-0.1*x(2)];

