Question: I've been working on a standard potential which I am trying to fit with a given model: ax2 - bx3 + lx4 The x and y values for the fit are generated from the code as well, the x values are generated by numpy.linspace. I have bounded the a,b,c parameters such that they are always positive. I needed the fit to mimic the data in such a way that at minimum, the height of the local maxima and the position of the global minima are accurate. Instead this is what I'm getting (blue is the actual data and the dashed line is the fitted one with the given model):  Here is the relevant part of my code: import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit from scipy.interpolate import * v =246. x = np.linspace(0.1, 270.1, 27001) Xdata = x/v def func(x,a,b,l): return (a*(x**2)) - (b*(x**3)) + ((l)*(x**4)) temp = np.linspace(80.,80.,1) b = np.zeros_like(temp) a = np.zeros_like(temp) l = np.zeros_like(temp) Vfit = pot_giver(temp[0]) #External func. Ydata = (Vfit - Vfit[0])/(pow(v,4)) popt, pcov = curve_fit(func, Xdata, Ydata,bounds=((0.,0.,0.), (np.inf,np.inf,np.inf))) b[0],a[0],l[0] = popt plt.plot(Xdata,Ydata) plt.plot(Xdata, func(Xdata, *popt), 'g--') plt.show()  I am doing this as an array because I am varying the temp over a wide range, this singular element 'temp' was created for troubleshooting. The Xdata and Ydata are here on the first and second columns respectively: Data

Answer: You can also use scipy.optimize.leastsq Here, you can also impose constraints (through a penalty added to the residuals when they are not satisfied). You asked for the following constraints:  the maximum of the function is correct;  the position of the global minimum is correct.   It's a little sensitive to how you impose the constraints, but this just about gets those. import math import numpy as np import matplotlib.pyplot as plt from scipy.optimize import leastsq Xdata, Ydata = np.loadtxt( "Pot_data_online.txt", unpack=True ) ymax = Ydata.max() ymin = Ydata.min() xymin = Xdata[0] for x, y in zip( Xdata, Ydata ): if y == ymin: xymin = x # Fitting function def f( x, a, b, l ): return a * x ** 2 - b * x ** 3 + l * x ** 4 # Residuals, including a penalty to get the constraints def residuals( p, x, y ): a, b, l = p x1 = ( 3 * b - math.sqrt( max( 9 * b ** 2 - 32 * a * l, 0 ) ) ) / ( 8 * l ) x2 = ( 3 * b + math.sqrt( max( 9 * b ** 2 - 32 * a * l, 0 ) ) ) / ( 8 * l ) error = y - f( x, a, b, l ) penalty = 1e6 * ( abs( f(x1,a,b,l) / ymax - 1 ) + abs( x2 / xymin - 1 ) ) return abs( error ) + penalty popt, pcov = leastsq( func=residuals, x0=(0.01,0.3,0.01), args=(Xdata,Ydata) ) a, b, l = popt x1 = ( 3 * b - math.sqrt( max( 9 * b ** 2 - 32 * a * l, 0 ) ) ) / ( 8 * l ) x2 = ( 3 * b + math.sqrt( max( 9 * b ** 2 - 32 * a * l, 0 ) ) ) / ( 8 * l ) print( "a, b, l = ", a, b, l ) print( "Constraint on maximum: ymax , fmax = ", ymax , f( x1, a, b, l ) ) print( "Constraint on position of minimum: xymin, xfmin = ", xymin, x2 ) plt.plot( Xdata,Ydata, 'b-' ) plt.plot( Xdata, f( Xdata, *popt ), 'g--') plt.show()  Output: a, b, l = 0.026786218092281242 0.0754627163616853 0.04481075395203686 Constraint on maximum: ymax , fmax = 0.0007404005868381408 0.0007404009309691355 Constraint on position of minimum: xymin, xfmin = 0.947308 0.947621778120905