program fskan
	implicit real *8 (a-h,o-z)
	real *8 eta(0:2000),f(0:2000),df(0:2000),ddf(0:2000)
	real *8 m,deta
	real *8 a2(500),dfl(500)
	write (6,1)
1	format('Input m')
	read(5,*)m
	b=2*m/(m+1)
	write(6,2)b
2	format(' b = ',f12.5)
	dflim=1.
c
c	search in range -0.1 < a2 <0.9 for value 
c	which has df->dflim as eta->infinity
c
	write(6,3)
3	format('Input search range:a2start, a2end')
	read(5,*)a2start,a2end
	na2=101
	da2=(a2end-a2start)/(na2-1)
	do k=1,na2
		a2(k)=a2start+(k-1)*da2
		deta=0.01
		n=1000
		call rkfskan(deta,n,eta,f,df,ddf,b,a2(k))
		dfl(k)=df(n)
		write(6,10)a2(k),dfl(k),ddf(n)
10	format(3f15.6)
	enddo
	write(8,12)(a2(k),dfl(k),k=1,na2)
12	format(f15.7,e16.7)
c
c	search for dflim amongst dfl values
c
	limfound=0
	do k=2,na2
		if((dfl(k)-dflim)*(dflim-dfl(k-1)).ge.0.)then
			limfound=1
			fac=(dflim-dfl(k-1))/(dfl(k)-dfl(k-1))
			a2found=(1-fac)*a2(k-1)+fac*a2(k)
			write(6,4)a2found
4			format('a2 value = ',f15.6)
6	format('Limit calculated = ',f15.6, 'limit sought = ',f15.6)
		endif
	enddo
	if(limfound.eq.1)then
	                call rkfskan(deta,n,eta,f,df,ddf,b,a2found)		
		        write(7,5)(eta(i),df(i),i=0,n)
5		format('data',/,(2f15.6))
		write(6,6)df(n),dflim
	else
		write(6,7)
7	format('No a2 value found')
	endif
	stop
	end
	subroutine rkfskan(deta,n,eta,f,df,ddf,b,a2)
	implicit real *8 (a-h,o-z)
	real *8 eta(0:2000),f(0:2000),df(0:2000),ddf(0:2000),a2
		eta(0)=0.
		f(0)=0
		df(0)=0
		ddf(0)=a2
		do i=1,n
		  fi=f(i-1)
		  dfi=df(i-1)
	  	  ddfi=ddf(i-1)

		  fa=deta*dfi
		  dfa=deta*ddfi
		  ddfa=deta*(b*(dfi*dfi-1.)-fi*ddfi)

		  fi=f(i-1)+fa/2
		  dfi=df(i-1)+dfa/2
		  ddfi=ddf(i-1)+ddfa/2

		  fb=deta*dfi
		  dfb=deta*ddfi
		  ddfb=deta*(b*(dfi*dfi-1.)-fi*ddfi)

		  fi=f(i-1)+fb/2
		  dfi=df(i-1)+dfb/2
		  ddfi=ddf(i-1)+ddfb/2

		  fc=deta*dfi
		  dfc=deta*ddfi
		  ddfc=deta*(b*(dfi*dfi-1.)-fi*ddfi)

		  fi=f(i-1)+fc
		  dfi=df(i-1)+dfc
		  ddfi=ddf(i-1)+ddfc

		  fd=deta*dfi
		  dfd=deta*ddfi
		  ddfd=deta*(b*(dfi*dfi-1.)-fi*ddfi)
		  eta(i)=eta(i-1)+deta
		  f(i)  =  f(i-1)+( fa+2* fb+2* fc+ fd)/6
		  df(i) = df(i-1)+(dfa+2*dfb+2*dfc+dfd)/6
		  ddf(i)=ddf(i-1)+(ddfa+2*ddfb+2*ddfc+ddfd)/6
		enddo
	return
	end