module ode
implicit none
contains

function rkstep12(f,t,y,h,err) result (yh)
real*8::t,y(:),h,err,yh(size(y))
interface
	function f(t,y) result (df)
	real*8::t,y(:),df(size(y))
	end function
end interface
real*8::k0(size(y)),k1(size(y))
k0=f(t,y)
k1=f(t+h/2,y+k0*h/2)
yh=y+k1*h
err=norm2(k1-k0)*h/2
end function rkstep12

function driver(f,a,y,b,acc,eps) result (path)
real*8::a,y(:),b,eps,acc
interface
	function f(t,y) result (df)
	real*8::t,y(:),df(size(y))
	end function
end interface
real*8,allocatable::path(:,:),tmp(:,:),yh(:)
real*8::t,h,tol,err
integer::m,n=1,k=100
m=size(y)
allocate(path(k,m+1))
t=a;h=(b-a)/20
path(n,1)=t
path(n,2:)=y
do while(n<=k)
	if(t>=b) then; exit ; end if
	if(t+h>b)then; h=b-t; end if
	yh=rkstep12(f,t,y,h,err)
	tol=(acc+norm2(yh)*eps)*sqrt(h/(b-a))
	if(err<tol)then
		n=n+1
		t=t+h
		y=yh
		if (n.eq.k) then
			tmp=path
			deallocate(path)
			k=k*2
			allocate(path(k,m+1))
			path(1:n,:)=tmp(:,:)
		end if
		path(n,1)=t
		path(n,2:)=y(:)
	end if
	if(err.eq.0)then; h=h*2
	else;             h=h*(tol/err)**0.25d0*0.95d0
	end if
end do
tmp=path
deallocate(path)
allocate(path(n,m+1))
path(:,:)=tmp(1:n,:)
end function driver

end module ode