How do I bootstrap a fitting function?

I have a fitting function below which fits a lorentzian to my data and I essentially want to firstly return the amplitude of my fit (perhaps more reliably than I have) and then secondly perform a bootstrap in order to calculate the standard error on my amplitude. However I have tried some bootstrap routines but either the value returned is way off what I expect or it isn’t able to fit anything.

This is my fitting function/routine so far :

def lorentzian(x, x0, a, gam ):
    return a * gam**2 / ( gam**2 + ( x - x0 )**2)

def multi_lorentz( x, params ):
    off = params[0]
    paramsRest = params[1:]
    assert not ( len( paramsRest ) % 3 )
    return off + sum( [ lorentzian( x, *paramsRest[ i : i+3 ] ) for i in range( 0, len( paramsRest ), 3 ) ] )

def res_multi_lorentz( params, xData, yData ):
    diff = [ multi_lorentz( x, params ) - y for x, y in zip( xData, yData ) ]
    return diff


x_data , y_data , name,y_data_2,bmi_err = get_y(combined_data,combined_info,36033)
x_peak = -26
x_range = [x_peak-25, x_peak+25]

# Filter data within the specified x_range
x_range_indices = np.where((x_data >= x_range[0]) & (x_data <= x_range[1]))[0]

if len(x_range_indices) == 0:
    raise ValueError("No data points within the specified x_range.")

x_range_data = x_data[x_range_indices]
y_range_data = y_data_2[x_range_indices]

bmi_error_range = bmi_err[x_range_indices]
x_fit = np.linspace(min(x_range_data),max(x_range_data),300)

xData, yData = x_range_data,y_range_data

generalWidth = 45

yDataLoc = yData
startValues = [ max( yData ) ]
counter = 0

while max( yDataLoc ) - min( yDataLoc ) > .1:
    counter += 1
    if counter > 20: ### max 20 peak...emergency break to avoid infinite loop
        break
    minP = np.argmin( yDataLoc )
    minY = yData[ minP ]
    x0 = -26
    startValues += [ x0, minY - max( yDataLoc ), generalWidth ]
    
    def model(x,*params):
        return multi_lorentz(x,params)
    
    
    popt,_ = curve_fit(model,xData,yData,p0 = startValues,maxfev= 500000,sigma = bmi_error_range )
    yDataLoc = [ y - multi_lorentz( x, popt ) for x,y in zip( xData, yData ) ]



testData = [ multi_lorentz(x, popt ) for x in x_fit ]

min_x_fit = x_fit[np.argmin(testData)]
max_x_fit = x_fit[np.argmax(testData)]
amplitude_at_min = multi_lorentz(min_x_fit,popt)
amplitude_at_max = multi_lorentz(max_x_fit,popt)
AMP = amplitude_at_min - amplitude_at_max