function harmonic;

x0=[10.;0.];
options=odeset('RelTol',1.e-6,'AbsTol',1.e-6)
[t,x]=ode45(@rpp,[0 100], x0,options);
ep=x(:,1).^2
ek=x(:,2).^2
etot=ep+ek
plot (t,ep(:),t,ek(:),t,etot(:),'-g')
function f=rpp(t,x)
f=[x(2);-x(1)]
