c optimization see Elementary numerical analysis,Conte& de Boor p.211
c      f(x,y)= x1**2 +x2**2
      f(x1,x2)=x1**3 + x2**3-2.*x1**2+3.*x2**2-8.
      dfdx1(x1,x2)=(f(x1+epsi,x2)-f(x1,x2))/epsi
      dfdx2(x1,x2)=(f(x1,x2+epsi)-f(x1,x2))/epsi
      dfdt(t)=(F(x1-(t+epsi)*dfdx1(x1,x2),x2-(t+epsi)*dfdx2(x1,x2) )-
     $ F(x1-t*dfdx1(x1,x2),x2-t*dfdx2(x1,x2)))/epsi
      dfdt2(t)=(dfdt(t+epsi)-dfdt(t) )/epsi
      x1=1.
      x2=-1.
      epsi=1.e-3
      print 100,x1,x2,f(x1,x2)
      print*,'dfdx1,dfdx2=',dfdx1(x1,x2),dfdx2(x1,x2)
      do 10 i=1,2
      t=1.
c the roots of the first derivative is calculated next it involves both
c  first derivative dfdt(t)  and the seciond dfdt2(t)
      do 20 j=1,10
      t=t-dfdt(t)/dfdt2(t)
20    continue
      z1=x1-t*dfdx1(x1,x2)
      z2=x2-t*dfdx2(x1,x2)
       print*,'z1,z2, f(z1,z2)=',z1,z2,f(z1,z2)
      x1=z1-t*dfdx1(z1,z2)
      x2=z2-t*dfdx2(z1,z2)
10    continue
100   format(1x,'x1,x2,f(x1,x2)=',3(3x,e10.3))
      stop
      end
      
      RUN
       x1,x2,f(x1,x2)=    0.100E+01   -0.100E+01   -0.700E+01
 dfdx1,dfdx2= -0.998973787 -2.99978232
 z1,z2, f(z1,z2)=  1.33253956 -0.00142905989 -9.1851778
 z1,z2, f(z1,z2)=  1.33233035 -0.00577284815 -9.18508339