clear all
close all
y0 = [0.5,2,10];
tspan = [0,15];
f = @(t,y) [-2*y(2)*y(1);...
    y(1)^2+y(3)^2-y(2)^2-1;...
    -2*(y(1)+y(2))*y(3)];
opzioni.InitialStep = 0.001;
opzioni.AbsTol = 1e-2;
opzioni.RelTol = 1e-2;
[tout,yout] = rk23(f,tspan,y0,opzioni);
E=(yout(:,1).^2+yout(:,2).^2+...
    yout(:,3).^2+1)./(2*yout(:,1));
plot(tout,E,'.')