c   8.Newton法で連立方程式を解く
c     森川 靖大 , 22980138
c
      program ex8
c
      implicit none
      integer :: n,i,j
      real :: a(10),b(10),x0,y0,xk1,yk1,epsilon
      epsilon = 1.0e-4    ! 収束精度

      do i=1,2
       if (i==1) then

c-------第1式の係数
        data a/ 0.0e0, 0.0e0, 0.0e0, 0.0e0, 1.0e0,
     -          0.0e0, 1.0e0, 0.0e0, 0.0e0,-1.0e0/

c-------第2式の係数
        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

c-------第1式の係数
        a(1) = 0.0e0           !第1式の係数
        a(2) = 0.0e0
        a(3) = 0.0e0
        a(4) = 0.0e0
        a(5) = 1.0e0
        a(6) = 0.0e0
        a(7) = 1.0e0
        a(8) = 0.0e0
        a(9) = 0.0e0
        a(10)=-1.0e0
c        data a/ 0.0e0, 0.0e0, 0.0e0, 0.0e0, 1.0e0,
c     -          0.0e0, 1.0e0, 0.0e0, 0.0e0,-1.0e0/

c-------第2式の係数
        b(1) = 0.0e0           !第2式の係数
        b(2) = 0.0e0
        b(3) = 0.0e0
        b(4) = 0.0e0
        b(5) = 0.0e0
        b(6) = 0.0e0
        b(7) = 1.0e0
        b(8) = 1.0e0
        b(9) = 0.0e0
        b(10)= 0.0e0
c        data b/ 0.0e0, 0.0e0, 0.0e0, 0.0e0, 0.0e0,
c     -          0.0e0, 1.0e0, 1.0e0, 0.0e0, 0.0e0/

       end if
       do j=1,2
        if  (j==1)  then
         x0=-0.1e0           ! 初期値
         y0=-1.0e0
        else if (j==2) then
         x0= 0.1e0           ! 初期値
         y0= 1.0e0
        endif

        call newtongragh(a,b,x0,y0,xk1,yk1,n,epsilon,i,j)
       enddo
      enddo
      stop
      end

c-----Newton法で2つのx,yの3次方程式の
c-----解と反復回数を求めるサブルーチン
      subroutine newtongragh(a,b,x0,y0,xk1,yk1,n,epsilon,i,j)
       implicit none
       integer :: n,i,j
       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)

       if (i==1 .and. j==1)then
       write(10,*)xk,yk
       else if (i==1 .and. j==2)then
       write(11,*)xk,yk
       else if (i==2 .and. j==1)then
       write(12,*)xk,yk
       else if (i==2 .and. j==2)then
       write(13,*)xk,yk
       end if

       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

       if (i==1 .and. j==1)then
       write(10,*)xk,yk
       else if (i==1 .and. j==2)then
       write(11,*)xk,yk
       else if (i==2 .and. j==1)then
       write(12,*)xk,yk
       else if (i==2 .and. j==2)then
       write(13,*)xk,yk
       end if

          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

      if (i==1 .and. j==1)then
      write(10,*)xk1,yk1
      else if (i==1 .and. j==2)then
      write(11,*)xk1,yk1
      else if (i==2 .and. j==1)then
      write(12,*)xk1,yk1
      else if (i==2 .and. j==2)then
      write(13,*)xk1,yk1
      endif

      return
      end

c-----関数 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-----関数 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-----関数 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-----関数 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-----関数 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-----関数 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