c   8.Newton$BK!$GO"N)J}Dx<0$r2r$/(B
c     $B?9@n(B $BLwBg(B , 22980138
c
      program ex8
c
      implicit none
      integer :: n,i,j
      real :: a(10),b(10),c(10),d(10),x0,y0,xk1,yk1,epsilon
      epsilon = 1.0e-4    ! $B<}B+@:EY(B

      write(6,*)'Newton$BK!$GO"N)J}Dx<0$r2r$/(B'
      write(6,300) epsilon
      do i=1,2
       if (i==1) then
        write(6,*)''
        write(6,*)'(1)'
        write(6,*)'$BO"N)J}Dx<0(B x**2+y**2-1=0'
        write(6,*)'           x**3-y=0     $B$N2r(B'

c-------$BBh(B1$B<0$N78?t(B
        data a/ 0.0e0, 0.0e0, 0.0e0, 0.0e0, 1.0e0,
     -          0.0e0, 1.0e0, 0.0e0, 0.0e0,-1.0e0/
       
c-------$BBh(B2$B<0$N78?t(B
        data b/ 1.0e0, 0.0e0, 0.0e0, 0.0e0, 0.0e0,
     -          0.0e0, 0.0e0, 0.0e0,-1.0e0, 0.0e0/
       
       else if (i==2) then
        write(6,*)''
        write(6,*)'(2)'
        write(6,*)'$BO"N)J}Dx<0(B x**2+y**2-1=0'
        write(6,*)'           x+y**2=0     $B$N2r(B'

c-------$BBh(B1$B<0$N78?t(B
        data c/ 0.0e0, 0.0e0, 0.0e0, 0.0e0, 1.0e0,
     -          0.0e0, 1.0e0, 0.0e0, 0.0e0,-1.0e0/
       
c-------$BBh(B2$B<0$N78?t(B
        data d/ 0.0e0, 0.0e0, 0.0e0, 0.0e0, 0.0e0,
     -          0.0e0, 1.0e0, 1.0e0, 0.0e0, 0.0e0/
       
       end if     
       do j=1,2
        if  (j==1)  then
         x0=-0.1e0           ! $B=i4|CM(B
         y0=-1.0e0
        else if (j==2) then
         x0= 0.1e0           ! $B=i4|CM(B
         y0= 1.0e0
        end if
        write(6,400) x0,y0
        
        if (i==1) then
          call newton(a,b,x0,y0,xk1,yk1,n,epsilon)
        elseif (i==2) then
          call newton(c,d,x0,y0,xk1,yk1,n,epsilon)
        end if
        write(6,200) n      
        write(6,100)'x=',xk1
        write(6,100)'y=',yk1
        write(6,*)''
       enddo
      enddo
      stop
  100 format(a3,e18.9)
  200 format(1x,'$BH?I|2s?t(B',i5) 
  300 format(1x,'$B<}B+@:EY(B',1x,e18.9)
  400 format(1x,'$B=i4|CM(B(x,y)=(',f4.1,',',f4.1,')') 
      end

c-----Newton$BK!$G(B2$B$D$N(Bx,y$B$N(B3$B<!J}Dx<0$N(B
c-----$B2r$HH?I|2s?t$r5a$a$k%5%V%k!<%A%s(B
      subroutine newton(a,b,x0,y0,xk1,yk1,n,epsilon)
       implicit none
       integer :: n
       real :: f,g,dfx,dgx,dfy,dgy,fs,gs,dfxs,dgxs,dfys,dgys,
     -         a(10),b(10), x0,y0,xk,yk,xk1,yk1,epsilon
      
       n=1
       xk=x0
       yk=y0
       fs  =  f(a,xk,yk)
       gs  =  g(b,xk,yk) 
       dfxs=dfx(a,xk,yk)
       dgxs=dgx(b,xk,yk)
       dfys=dfy(a,xk,yk)
       dgys=dgy(b,xk,yk)

       xk1=xk-( fs*dgys-gs*dfys)/(dfxs*dgys-dfys*dgxs)
       yk1=yk-(-fs*dgxs+gs*dfxs)/(dfxs*dgys-dfys*dgxs)
  
        do while (sqrt((xk1-xk)**2+(yk1-yk)**2) >= epsilon)
          n=n+1
          xk=xk1
          yk=yk1

          fs  =  f(a,xk,yk)
          gs  =  g(b,xk,yk) 
          dfxs=dfx(a,xk,yk)
          dgxs=dgx(b,xk,yk)
          dfys=dfy(a,xk,yk)
          dgys=dgy(b,xk,yk)

          xk1=xk-( fs*dgys-gs*dfys)/(dfxs*dgys-dfys*dgxs)
          yk1=yk-(-fs*dgxs+gs*dfxs)/(dfxs*dgys-dfys*dgxs)
        enddo

      return
      end

c-----$B4X?t(B f(x,y)
      function f(a,xk,yk)
      implicit none
      real :: f,xk,yk,a(10)
      
      f= a(1)*xk**3 +a(2)*xk**2*yk +a(3)*xk*yk**2 +a(4)*yk**3
     -  +a(5)*xk**2 +a(6)*xk*yk +a(7)*yk**2 +a(8)*xk +a(9)*yk +a(10)

      return
      end

c-----$B4X?t(B g(x,y)
      function g(b,xk,yk)
      implicit none
      real :: g,xk,yk,b(10)
      
      g= b(1)*xk**3 +b(2)*xk**2*yk +b(3)*xk*yk**2 +b(4)*yk**3
     -  +b(5)*xk**2 +b(6)*xk*yk +b(7)*yk**2 +b(8)*xk +b(9)*yk +b(10)

      return
      end

c-----$B4X?t(B df(x,y)/dx
      function dfx(a,xk,yk)
      implicit none
      real :: dfx,xk,yk,a(10)

      dfx= 3.0e0*a(1)*xk**2 +2.0e0*a(2)*xk*yk +a(3)*yk**2
     -    +a(4)*0.0e0 +2.0e0*a(5)*xk +a(6)*yk +a(7)*0.0e0
     -    +a(8) +a(9)*0.0e0 +a(10)*0.0e0

      return
      end

c-----$B4X?t(B dg(x,y)/dx
      function dgx(b,xk,yk)
      implicit none
      real :: dgx,xk,yk,b(10)

      dgx= 3.0e0*b(1)*xk**2 +2.0e0*b(2)*xk*yk +b(3)*yk**2
     -    +b(4)*0.0e0 +2.0e0*b(5)*xk +b(6)*yk +b(7)*0.0e0
     -    +b(8) +b(9)*0.0e0 +b(10)*0.0e0

      return
      end

c-----$B4X?t(B df(x,y)/dy
      function dfy(a,xk,yk)
      implicit none
      real :: dfy,xk,yk,a(10)
      
      dfy= a(1)*0.0e0 +a(2)*xk**2 +2.0e0*a(3)*xk*yk
     -    +3.0e0*a(4)*yk**2 +a(5)*0.0e0 +a(6)*xk +2.0e0*a(7)*yk
     -    +a(8)*0.0e0 +a(9) +a(10)*0.0e0

      return
      end

c-----$B4X?t(B dg(x,y)/dy
      function dgy(b,xk,yk)
      implicit none
      real :: dgy,xk,yk,b(10)
      
      dgy= b(1)*0.0e0 +b(2)*xk**2 +2.0e0*b(3)*xk*yk
     -    +3.0e0*b(4)*yk**2 +b(5)*0.0e0 +b(6)*xk +2.0e0*b(7)*yk
     -    +b(8)*0.0e0 +b(9) +b(10)*0.0e0

      return
      end
