def Boettger(ax, layers, energies, area, tag):
    """
    Function that produces the surface energy result from the Boettger Method. 
    Call function for each relevant surface. 
    
    Arguments
    ---------
    ax: matplotlib.axes
        Figure axes to add the linear line to. 
    layers: iterable
        Sorted iterable of integers corresponding to the number of layers that 
        will make up the x-axis. 
    energies: iterable
        Sorted w.r.t. the number of layers and will make up the y-axis. 
    area: float
        Surface area of the relevant surface. 
    tag: str
        String indicating the method. Must be one of: "pbe", "ts", "mbd"
    
    """
    x, y = [], []
    c = {"pbe":"orange", "ts":"g", "mbd":"r"}
    for i in range(2, len(layers)):
        bulk = (energies[i] - energies[i-1]) / (layers[i] - layers[i-1])
        x.append(layers[i])
        y.append(16000*(energies[i] - layers[i] * bulk)/(2*area))
    ax.plot(x, y, c=c[tag], label="{} Boettger".format(
        tag.upper()), linewidth=7.0)
    return x, y