clear all
close all
tspan = [0,15];
y0 = [0.5;2;10];
tol = 1e-4;
options.RelTol = tol;
options.AbsTol = tol;
options.InitialStep = ((1+norm(y0))*tol)^(1/3);
[tout,yout] = rk23(@schrof,tspan,y0,options);
E = (yout(1,:).^2+yout(2,:).^2+yout(3,:).^2+1)./(2*yout(1,:));
plot(tout,E,'.')
xlabel('t')
ylabel('E(t)')
tspan = [0,1];
tol = 1e-10
options.RelTol = tol;
options.AbsTol = tol;
options.InitialStep = 1e-3;
[tout,youtrif] = rk23(@schrof,tspan,y0,options);
tolrange = [1e-5,1e-6,1e-7,1e-8];
counter = 0;
for tol = tolrange
  tol
  counter = counter+1;
  options.RelTol = tol;
  options.AbsTol = tol;
  options.InitialStep = ((1+norm(y0))*tol)^(1/3);
  [tout,yout] = rk23(@schrof,tspan,y0,options);
  err(counter) = norm(yout(:,end)-youtrif(:,end));
  timesteps(counter) = length(tout);
end
figure
loglog(timesteps,err,'*',timesteps,err(end)*(timesteps/timesteps(end)).^(-3))
legend('error','order 3')
xlabel('timesteps')
ylabel('error')