eqlbrm.f

C     astar & bstar: guesses in, optimized values out
      subroutine eqlbrm(biga, bigb, astar, bstar, n)
      external biga, bigb
        function f(astar)
          function g(a)
            function h(b)
            h = bigb(astar, b)
            end
          bstar = argmax(h, bstar, n)
          g = biga(a, bstar)
          end
        f = argmax(g, astar, n) - astar
        end
      astar = root(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

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

      function argmax(f, x0, n)
        function fprime(x)
        fprime = deriv1(f, x)
        end
      argmax = root(fprime, 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