Scipy.optimzie can't get the right fit for data with errors (polynomial function 3rd order)

Question

so I have a problem with fitting a polynomial function 3rd order to my measured data points. The fit I get looks good and finds the minimum quite well, but the parameters and their errors are very big and useless. I know aswell, that for example numpy.polynomial or Fityk (program for fitting) get the same results, so it's not the fault of scipy.optimize.

Here's a minimalistic example:

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

x = [113.5,129.9,134.9,140.0,143.4,145.0,146.0,146.5,147.0,147.4,148.0,148.4,148.9,149.4,149.9,150.4,151.0,151.1,151.4,151.5,151.7,152.0,152.2,152.4,152.6,152.8,152.9,153.1,153.4,153.5,153.8,154.0,154.2,154.4,154.6,154.7,155.0,155.2,155.4,155.5,155.8,156.0,156.2,156.3,156.5,156.7,156.9,157.2,157.3,157.5,157.8,157.9,158.1,158.4,158.6,158.8,159.0,159.2,159.4,159.6,159.8,160.0,160.4,160.9,161.4,161.9,162.5,162.9,163.4,163.9,164.4,164.9,165.4]
y = [9.924496174702561,5.1614753220580445,4.273935632622366,3.094738666201405,2.083925837022858,2.115431896180724,1.9113528763921894,1.879355666362959,1.679016362042712,1.7271227281276254,1.6609811891540538,1.5249993702844755,1.5129258663325422,1.500042284124191,1.3932211464841808,1.3886628878493803,1.2929918423598743,1.2874924482437715,1.2948585494988538,1.3075898522476532,1.2201992422487522,1.264328864187696,1.215889388646001,1.2319937419934728,1.1906273595713857,1.149716935155421,1.1723839505427636,1.1036830323762243,1.164896884175834,1.154067184533763,1.157669109948195,1.098427466779971,1.0951354571029084,1.141132347068189,1.1010769025947835,1.1737462215496937,1.2075760000300582,1.1521169161001228,1.145815193564697,1.208010688484708,1.1806288211421971,1.2094270267668206,1.1926204409804952,1.2371404098805436,1.2505603552589275,1.2597590363411229,1.2340000589814226,1.245960988062385,1.2813709291729314,1.2858127915443631,1.3309207448989704,1.3414497752361316,1.3366393982318492,1.3769180494934625,1.374487109711036,1.4190515887645334,1.4087845418301346,1.4670046473870004,1.5009920167035113,1.516433575597311,1.5581877808935867,1.5639034017453048,1.6432985329201213,1.6806702927177568,1.7603402605713827,1.7726684359607014,1.8425660780655149,1.9397009344863552,2.034447328713384,2.0692854680277253,2.13766242918302,2.2265883882041977,2.2973788056787914]
dy = [0.06511217866171501,0.018759889812228713,0.014960007743290624,0.008436799186629396,0.004641280772920625,0.003806943299498247,0.0033080468027512377,0.0031292816259633075,0.0034897125415282813,0.0033805656175818323,0.0032118341411569084,0.003456776531758499,0.0035130943741483035,0.003844669995096197,0.003705111855426419,0.004265911409700588,0.004119463555766376,0.004333160800489313,0.004675586935753023,0.004755844981466414,0.004626063907657897,0.00484418171028472,0.004975598660818002,0.004794722872056576,0.004935523410534309,0.0047409725117238295,0.005285678816054232,0.005647757266973509,0.005233111822412779,0.005077410240144215,0.0055954921053261355,0.006005301592808242,0.005891498159940475,0.005912222990581101,0.006222905717368581,0.00589746079266117,0.006011573684638575,0.005788227543880991,0.005962188009592335,0.006225357597091497,0.006222985355473735,0.006541769810610297,0.006870025948951895,0.006653608687637577,0.00639622489493887,0.006369179745681897,0.006499053395970663,0.006461841584749265,0.006951496563160294,0.006868722120219894,0.007035010105556161,0.0068105211501542865,0.006300523665919444,0.006429542410031948,0.006718691462157191,0.00705522218620391,0.006628996527755684,0.007039161760455189,0.006971644656870013,0.007067453435110972,0.007052037797103016,0.006894427520333761,0.0068549023085317405,0.006951388538832804,0.007037404310773895,0.0068284487430299665,0.006551746301744071,0.0065459881397853656,0.0066975657448055306,0.006432309692561562,0.006310794127058911,0.006086817068652075,0.006067627534027548]
dx = 0.1

def fit(x,a,b,c,d):
    return a*x**3+b*x**2+c*x+d
    
params, cov = optimize.curve_fit(fit,x,y)
errors = np.sqrt(np.diag(cov))

print(params,"\n",errors)    

plt.errorbar(x,y,xerr=dx,yerr=dy,fmt='r.',label="datapoints")

xfit = np.arange(min(x),max(x),0.01)    #so the fit looks nice
plt.plot(xfit,fit(xfit,*params),label="optimize-fit")

plt.legend()
plt.show()

I'm sorry about the way of giving you the data points, but I didn't know how to give them to you another way.
The output of the code is:
[ 1.06504223e-04 -3.98864448e-02 4.69539864e+00 -1.64923300e+02] [4.69911453e-06 2.00337778e-03 2.82662482e-01 1.31848506e+01]
And a closed up part of the plot looks like this: plot by the given code.

I think the problem is, that I don't give optimize the errors of my variables, but if I add sigma with and without absolute_sigma it doesn't get better. Does anyone knows how to improve my code/fit? Maybe there is a better function I haven't thought of that suits the data better.

Btw this is a copied question from Stack Overflow (https://stackoverflow.com/questions/62360580/scipy-optimzie-cant-get-the-right-fit-for-data-with-errors-polynomial-function), but i thought since nobody answered there maybe someone here knows an answer.