

global omega;
k=10.;m=1.;
omega=k/m;
x0=[1.;0.];
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',1.);
[t,x]=ode45(@rpp,[0 4], x0,options);
epot=0.5*k*x(:,1).^2;ekin=0.5*m*x(:,2).^2;
ekokku=ekin+epot;
plot(t,[ekokku,epot,ekin]);
