c Numerical solution to the Dirac equation where
c   G=large component, F=small component
c  see also www.oocities.org/serienumerica
c by Reinaldo Baretti Machin
c Depto de Fisica UPR-Humacao
c reibaretti2004@yahoo.com
c www.oocities.org/reibaretti2004


      dimension F(0:10000), G(0:10000)
      real kapa
      v(x)=-z/x
      eminus(x)=erel-v(x)-c**2
      eplus(x)=erel-v(x)+c**2
      c=137.11
      z=30.
      kapa=-1.
      ei=-1.1*z**2/(2.)
      ef=-.9*z**2/(2.)
      ie=20
      de=(ef-ei)/float(ie)
      e=ei
      nx=2000
      print*,'edirac=',c**2*(sqrt(1.-(z/c)**2)-1.)
      print*,' '
      do 40 it=1,ie
      xf=-3.8*z/e
      dx=xf/float(nx)
      erel=e+c**2
      f(0)=0.
      g(0)=0.
      g(1)=1.
      f(1)=(Z/(2.*c))*g(1)
      do 10 i=2,nx
      x=dx*float(i)
      g(i)=g(i-1)+dx*(-kapa*g(i-1)/(x-dx)+(1./c)*eplus(x-dx)*f(i-1))
      f(i)=f(i-1)+dx*(kapa*f(i-1)/(x-dx)-(1./c)*eminus(x-dx)*g(i-1))
10    continue
      print 100,e,g(nx),f(nx)
      e=e+de
40    continue
100   format(1x,'e,g(large),f(small)=',3(3x,e11.4))
      stop
      end
      
      The asymptotic wavefunction changes sign when the the eigenvalue is       reached. 
       edirac= -455.519226

 e,g(large),f(small)=   -0.4950E+03    0.6762E+03    0.5459E+02
 e,g(large),f(small)=   -0.4905E+03    0.6099E+03    0.4904E+02
 e,g(large),f(small)=   -0.4860E+03    0.5411E+03    0.4333E+02
 e,g(large),f(small)=   -0.4815E+03    0.4697E+03    0.3746E+02
 e,g(large),f(small)=   -0.4770E+03    0.3956E+03    0.3140E+02
 e,g(large),f(small)=   -0.4725E+03    0.3187E+03    0.2517E+02
 e,g(large),f(small)=   -0.4680E+03    0.2389E+03    0.1875E+02
 e,g(large),f(small)=   -0.4635E+03    0.1560E+03    0.1213E+02
 e,g(large),f(small)=   -0.4590E+03    0.6989E+02    0.5314E+01
 e,g(large),f(small)=   -0.4545E+03   -0.1955E+02   -0.1714E+01
 e,g(large),f(small)=   -0.4500E+03   -0.1125E+03   -0.8962E+01
 e,g(large),f(small)=   -0.4455E+03   -0.2091E+03   -0.1644E+02
 e,g(large),f(small)=   -0.4410E+03   -0.3096E+03   -0.2415E+02
 e,g(large),f(small)=   -0.4365E+03   -0.4140E+03   -0.3210E+02
 e,g(large),f(small)=   -0.4320E+03   -0.5227E+03   -0.4032E+02
 e,g(large),f(small)=   -0.4275E+03   -0.6357E+03   -0.4880E+02
 e,g(large),f(small)=   -0.4230E+03   -0.7534E+03   -0.5755E+02
 e,g(large),f(small)=   -0.4185E+03   -0.8759E+03   -0.6659E+02
 e,g(large),f(small)=   -0.4140E+03   -0.1003E+04   -0.7594E+02
 e,g(large),f(small)=   -0.4095E+03   -0.1136E+04   -0.8559E+02