I'm trying to fit and plot a sigmoid curve fitted to my data. Here is my code:

import numpy as np
import pyplot
from scipy.optimize import curve_fit

def sigmoid(x, a, b):
    return 1.0 / (1.0   np.exp(-a*(x-b)))

xdata = [ 26457.28874429,  32614.0743368 ,  36375.76229015,  39371.21198926,
    41892.69310906,  44162.12906554,  46223.90776821,  48197.46060769,
    50022.99546123,  51848.03215704,  53568.91252972,  55354.78073918,
    57040.61643994,  58818.59692639,  60676.00000216,  62510.8110985 ,
    64586.59777512,  66562.62315247,  68742.97216949,  71078.33629888,
    73648.94810531,  76587.27871678,  79749.19546409,  83334.78854123,
    87282.81262044,  91925.69906738,  97551.59686331, 104496.29453508,
   112988.76665607, 124907.20130917, 144733.52923071, 193062.81222137]
ydata = [0.        , 0.57224724, 0.67521945, 0.76716544, 0.7535795 ,
   0.69637917, 0.81766149, 0.80906626, 0.71942446, 0.7670196 ,
   0.77549438, 0.80370093, 0.77160494, 0.78232104, 0.87796313,
   0.89326037, 1.05848115, 1.08499096, 1.13992674, 1.16009281,
   1.36193537, 1.58227848, 1.78731153, 1.98694295, 2.19359185,
   2.51098556, 2.92349108, 3.47826087, 3.82897412, 4.23472474,
   4.19622344, 5.19657584]

         
popt, pcov = curve_fit(sigmoid, xdata, ydata, method='dogbox', bounds=([min(ydata), min(xdata)],[max(ydata), max(xdata)]))
print(popt)

x = np.linspace(np.min(xdata), np.max(xdata))
y = sigmoid(x, *popt)

pyplot.plot(xdata, ydata, 'o', label='data')
pyplot.plot(x,y, label='fit')
pyplot.legend(loc='best')
pyplot.show()

However, I get an unsmoothed-squared-like curve like this:

I'm following this