real y(6),work(39),b(3),e(3),m
integer iwork(5)
common qm,b,e
external f
namelist /input/q,m,nt,dt,abserr,relerr,b,e,y
open(10,file="lorentz.in")
open(20,file="lorentz.dat")

read(10,input)
print input
qm=q/m
neqn=6
t=0
iflag=1

do i=1,nt
tout=t+dt

call rkf45(f,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
if(iflag.ne.2) then
iflag=2
endif

write(20,*) tout,y
enddo

stop
end
	
subroutine f(t,y,dy)
real y(6),dy(6),b(3),e(3)
common qm,b,e
!b(3)=(0.1*y(3))**4
dy(1)=y(4)
dy(2)=y(5)
dy(3)=y(6)
dy(4)=qm*(y(5)*b(3)-y(6)*b(2)+e(1))
dy(5)=qm*(y(6)*b(1)-y(4)*b(3)+e(2))
dy(6)=qm*(y(4)*b(2)-y(5)*b(1)+e(3))

return
end