function [t,y,H]=ode12(f,tspan,y0,tol);
% ode12 solves system of ODEs using adaptive integration
%   [t,y]=ode12(f,tspan,y0,tol) solves dy/dt=f(t,y) with initial
%   value y0 on the time interval tspan using adaptive integration
%   with Euler and Runge's method and error tolerance tol.
%   For didactical reasons we store the used stepsizes in H. 

dt=(tspan(2)-tspan(1))/10; H=dt;
y(:,1)=y0;
t(1)=tspan(1);
as=0.8; amin=0.2; amax=5;              % heuristic parameters
i=1;
while t(i)<tspan(2),
  k1=f(t(i),y(:,i));                      
  k2=f(t(i)+dt/2,y(:,i)+dt/2*k1);
  yn=y(:,i)+dt*k2;                     % Runge step
  yt=y(:,i)+dt*k1;                     % Euler step
  toli=max([norm(yn,inf) norm(y(:,i),inf)]); 
  toli=tol*(1+toli);                   % abs. and relative tol
  errest=norm(yn-yt);
  if errest<toli,                      % accept the step
    i=i+1;
    y(:,i)=yn;
    t(i)=t(i-1)+dt;
  end;
  dtopt=as*(dt*sqrt(toli/errest));     % estimate the new time step
  if dtopt<amin*dt,
    dt=amin*dt;
  elseif dtopt>amax*dt,
    dt=amax*dt;
  else
    dt=dtopt;
  end;
  if t(i)+dt>tspan(2),
    dt=tspan(2)-t(i); 
  end;
  H=[H;dt];
end;
t=t';
y=y';