def bisection(f,a,b,N):
'''Approximate solution of f(x)=0 on interval [a,b] by the bisection method.

Parameters
----------
f : function
    The function for which we are trying to approximate a solution f(x)=0.
a,b : numbers
    The interval in which to search for a solution. The function returns
    None if f(a)*f(b) >= 0 since a solution is not guaranteed.
N : (positive) integer
    The number of iterations to implement.

Returns
-------
x_N : number
    The midpoint of the Nth interval computed by the bisection method. The
    initial interval [a_0,b_0] is given by [a,b]. If f(m_n) == 0 for some
    midpoint m_n = (a_n + b_n)/2, then the function returns this solution.
    If all signs of values f(a_n), f(b_n) and f(m_n) are the same at any
    iteration, the bisection method fails and return None.

Examples
--------
>>> f = lambda x: x**2 - x - 1
>>> bisection(f,1,2,25)
1.618033990263939
>>> f = lambda x: (2*x - 1)*(x - 3)
>>> bisection(f,0,1,10)
0.5
'''
if f(a)*f(b) >= 0:
    print("Bisection method fails.")
    return None
a_n = a
b_n = b
for n in range(1,N+1):
    m_n = (a_n + b_n)/2
    f_m_n = f(m_n)
    if f(a_n)*f_m_n < 0:
        a_n = a_n
        b_n = m_n
    elif f(b_n)*f_m_n < 0:
        a_n = m_n
        b_n = b_n
    elif f_m_n == 0:
        print("Found exact solution.")
        return m_n
    else:
        print("Bisection method fails.")
        return None
return (a_n + b_n)/2

import pandas as pd
   df = pd.read_excel("directory.file.xlsx")
   v= df
   v

   D    P   h   O   Q   SD  LT
   80   27  0.37    50  2   3   1.51

from scipy.stats import norm
     import numpy as np
     f = lambda x: norm.ppf(1- (v['Q']*v['P']*v['h'])/(2*v['D']*x))*
     v['SD']*np.sqrt(v['LT'])*v['P']*v['h']+
     np.sqrt(2*v['D']*v['O']*v['h']*v['P'])-
     v['D']*x*(((-(norm.ppf(1-(v['Q']*v['h']*
     v['P'])/(2*x*v['D']))))*(1-norm.cdf((norm.ppf(1-
     (v['Q']*v['h']*v['P'])/(2*x*v['D']))),loc=0,scale=1))
     +(norm.pdf((norm.ppf(1-
     (v['Q']*v['h']*v['P'])/(2*x*v['D']))),loc=0,scale=1)))
     *v['SD']*np.sqrt(v['LT'])-v['Q'])/v['Q']*-1
     b = 10
     a = 1
     N = 100
     v['f'] = v.apply(bisection(f,a,b,N),axis=0)
     v

from scipy.stats import norm
    import numpy as np
    f = lambda x: norm.ppf(1- (2*27*0.37)/(2*80*x))*3*np.sqrt(1.51)*
    27*0.37+np.sqrt(2*80*50*0.37*27)-80*x*
    (((-(norm.ppf(1-(2*0.37*27)/(2*x*80))))*
    (1-norm.cdf((norm.ppf(1-(2*0.37*27)/(2*x*80))),loc=0,scale=1))+
    (norm.pdf((norm.ppf(1-(2*0.37*27)/(2*x*80))),loc=0,scale=1)))*
    3*np.sqrt(1.51)-2)/2*-1
    result = bisection(f,1,10,100)
    print(result)

   Found exact solution.
   4.503647631997683

from scipy.stats import norm
     import numpy as np
     f = lambda x: norm.ppf(1- (v['Q']*v['P']*v['h'])/(2*v['D']*x))*
     v['SD']*np.sqrt(v['LT'])*v['P']*v['h']+
     np.sqrt(2*v['D']*v['O']*v['h']*v['P'])-
     v['D']*x*(((-(norm.ppf(1-(v['Q']*v['h']*
     v['P'])/(2*x*v['D']))))*(1-norm.cdf((norm.ppf(1-
     (v['Q']*v['h']*v['P'])/(2*x*v['D']))),loc=0,scale=1))
     +(norm.pdf((norm.ppf(1-
     (v['Q']*v['h']*v['P'])/(2*x*v['D']))),loc=0,scale=1)))
     *v['SD']*np.sqrt(v['LT'])-v['Q'])/v['Q']*-1
     b = 10
     a = 1
     N = 100
     v['f'] = v.apply(bisection(f,a,b,N),axis=0)
     v