rk4adapt(f, span, y, epsi)

function y = rk4adapt(f, span, y, epsi)

span = sort(span);
steps=abs(diff(span));
for i=1:length(span)-1
   y = [y  rk4interne(f, span(i), y(:,i), steps(i), epsi)];
end;

function y = rk4interne(f, x, y, h, epsi)
y1 = rk4(f, x, y, h);
y2 = rk4(f, x, y, h/2);
y2 = rk4(f, x+h/2, y2, h/2);
if abs(y1-y2)<15*epsi
   y = y2;
else
   y = rk4interne(f, x, y, h/2, epsi);
   y = rk4interne(f, x+h/2, y, h/2, epsi);
end;