Minimizing Chi-Squared in Matplotlib

Pretty much, I am trying to use Chi-Squared to find the best fit for my model. My model is dependent on two fixing parameters as seen in the code: p0 and r0. I need to find the values of these two parameters in which the error between the model and the observed data is minimal. While I know how to do this on a basic level with simply minimizing the value of a product of say the numbers 2 and 5, I don't know how to minimize with arrays of data. Here is my code (Ignore the bottom half of the code which is just plotting the graph):

from scipy.optimize import*
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import scipy.special as special
from scipy.integrate import quad
from matplotlib.lines import Line2D
import matplotlib.ticker as ticker

#number of measurements
am = 27



#Observed Data 
v_obs = [
#Some data

]

#Fixing Parameters

p0 = 2.95*10**8
r0 = 1.8

#integral
def integrand(r,p0,r0):
    return (p0 * r0**3)/((r+r0)*((r**2)+(r0**2)))*4*3.1415926535*r**2
integrals = []
for i in r:
     integrals.append(quad(integrand, 0 ,i,args=(p0,r0)))


#function/Model
functions = []


for x in range (0,am):
    k = integrals[x][0]
    i = r[x]
    functions.append(np.sqrt((4.302*10**(-6)*k)/i))


#Plot the function/model
plt.plot(r,functions,color='red')




#Gas velocity (Don't need to worry about this)
v_gas = (
#Some Data
)


#Error Bars for Obs (with errors varying for each measurement)
#I apologize for the messy code

errors = []
errors.append(3.62)
errors.append(4.31)
errors.append(3.11)
errors.append(5.5)
errors.append(3.9)
errors.append(3.5)
errors.append(2.7)
errors.append(2.5)
errors.append(2.3)
errors.append(2.1)
errors.append(2.3)
errors.append(2.6)
errors.append(3.1)
errors.append(3.2)
errors.append(3.2)
errors.append(3.1)
errors.append(2.9)
errors.append(2.8)
errors.append(3.3)
errors.append(3.9)
errors.append(4.4)
errors.append(4.8)
errors.append(4.4)
errors.append(3.35)
errors.append(3.95)
errors.append(4.44)
errors.append(5.36)

for x in range (0,am):
    plt.errorbar(r[x],v_obs[x],yerr = errors[x],capsize = 2, capthick = 2, linestyle = 'None', color = 'green')



#Velocity v.s Radius
plt.scatter(r, v_obs, s = 8)
plt.plot(r,v_gas, '--',r,v_disk,'--')

ax = plt.gca()
ax.text(.2, 65, 'ps ' + str(p0) +  '\nrs ' + str(r0), style='italic')
        #bbox={'facecolor':'red', 'alpha':0.5, 'pad':10})




#Legend the Plot
plt.title("NGC2976 Rotational Curve")
plt.xlabel("Radius (kpc)")
plt.ylabel("Velocity (km/s)")

#fig = plt.figure()                                                               
#ax = fig.add_subplot(1,1,1)                                                      

# major ticks every 20, minor ticks every 5                                      
majory_ticks = np.arange(0, 120, 10)                                              
minory_ticks = np.arange(0, 120, 2)

majorx_ticks = np.arange(0,3,.3)
minorx_ticks = np.arange(0,3,.06)

ax.set_yticks(majory_ticks)                                                       
ax.set_yticks(minory_ticks, minor=True)                                           
ax.set_xticks(majorx_ticks)                                                       
ax.set_xticks(minorx_ticks, minor=True)                                           

# and a corresponding grid                                                       
ax.grid(which='both')                                                            

# or if you want differnet settings for the grids:                               
ax.grid(which='minor', alpha=0.2)                                                
ax.grid(which='major', alpha=0.5)

plt.show()