eqlbrm7.f

      function eqlbrm_f(astar)
      external biga, bigb
        function g(a)
          function h(b)
          h = bigb(astar, b)
          end
        bstar = argmax(h, bstar, n)
        g = biga(a, bstar)
        end
      eqlbrm_f = argmax(g, astar, n) - astar
      end

C     astar & bstar: guesses in, optimized values out
      subroutine eqlbrm(biga, bigb, astar, bstar, n)
      astar = root(eqlbrm_f, astar, n)
      end

      function root(f, x0, n)
      x = x0
      do 1669 i = 1, n
      call deriv2(f, x, y, yprime)
 1669 x = x - y / yprime
      root = x
      end

      function root_1(f1, x0, n)
      x = x0
      do 1669 i = 1, n
      call deriv2_1(f1, x, y, yprime)
 1669 x = x - y / yprime
      root_1 = x
      end

      subroutine deriv2(f, x, y, yprime)
      external f
      adf(x)
      y = f(x)
      end adf(yprime = tangent(y))
      end

      subroutine deriv2_1(f1, x, y, yprime)
      adf(x)
      y = argmax_fprime(x, f1)
      end adf(yprime = tangent(y))
      end

      function argmax_fprime(x, f)
      argmax_fprime = deriv1(f, x)
      end

      function argmax(f, x0, n)
      argmax = root_1(f, x0, n)
      end

      function deriv1(f, x)
      external f
      adf(x)
      y = f(x)
      end adf(deriv1 = tangent(y))
      end

      function gmbiga(a, b)
      price = 20 - 0.1*a - 0.1*b
      costs = a*(10 - 0.05*a)
      gmbiga = a*price - costs
      end

      function gmbigb(a, b)
      price = 20 - 0.1*b - 0.0999*a
      costs = b*(10.005 - 0.05*b)
      gmbigb = b*price - costs
      end

      program main
      read *, astar
      read *, bstar
      read *, n
      call eqlbrm(gmbiga, gmbigb, astar, bstar, n)
      print *, astar, bstar
      end