real y(6),work(39),b(3),e(3)
	integer iwork(5)
	external func
	common emq,e,b,oq,om,beta
	namelist /input/ nt,dt,neqn,relerr,abserr,eq,em,oq,om,
     *  beta,e,b,y
	open (10,file="lorenz.in")
	read (10,input) 
        emq=eq*oq/(em*om)
        t=0
        tou=t
        iflag=1
        call rkf45(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
	do it=1,nt
	iflag=2
	tout=t+dt
        call rkf45(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
	t=tout
	print *,y(1:3)
	enddo
	stop
	end
	
	subroutine func(t,y,dy)
	common emq,e,b,oq,om,beta
	real y(6),dy(6),e(3),b(3)
	vv=sqrt(y(4)**2+y(5)**2+y(6)**2)
	dy(1)=y(4)
	dy(2)=y(5)
	dy(3)=y(6)
	dy(4)=(y(5)*b(3)-y(6)*b(2)+e(1))*emq-beta*vv*y(4)
	dy(5)=(y(6)*b(1)-y(4)*b(3)+e(2))*emq-beta*vv*y(5)
	dy(6)=(y(4)*b(2)-y(5)*b(1)+e(3))*emq-beta*vv*y(6)
	return
	end