program fskan
c
c	search for solutions setting a2(0) and varying b
	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 b(500),dfl(500)
	write (6,1)
1	format('Input a2')
	read(5,*)a2
	dflim=1.
c
c	search in range of beta for value 
c	which has df->dflim as eta->infinity
c
	write(6,3)
3	format('Input search range:b start, bend')
	read(5,*)bstart,bend
	nb=201
	db=(bend-bstart)/(nb-1)
	do k=1,nb
		b(k)=bstart+(k-1)*db
		deta=0.01
		n=2000
		call rkfskan(deta,n,eta,f,df,ddf,b(k),a2)
		dfl(k)=df(n)
		write(6,10)b(k),dfl(k),ddf(n)
10	format(3f15.6)
	enddo
	write(8,12)(b(k),dfl(k),k=1,nb)
12	format(f15.7,e16.7)
c
c	search for dflim amongst dfl values
c
	limfound=0
	do k=2,nb
		if((dfl(k)-dflim)*(dflim-dfl(k-1)).ge.0.)then
			limfound=1
			fac=(dflim-dfl(k-1))/(dfl(k)-dfl(k-1))
			bfound=(1-fac)*b(k-1)+fac*b(k)
			write(6,4)bfound
4			format('beta value = ',f15.6)
	                call rkfskan(deta,n,eta,f,df,ddf,bfound,a2)		
		        write(7,5)(eta(i),df(i),i=0,n)
5		format('data',/,(2f15.6))
		write(6,6)df(n),dflim
6	format('Limit calculated = ',f15.6, 'limit sought = ',f15.6)
		endif
	enddo
	if(limfound.eq.1)then
	else
		write(6,7)
7	format('No b 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