How to fit and plot an smoothed sigmoid function?

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 answer

So there were two problems I encountered:

1. Your function was inappropriate to your data. Note your function would never be greater than 1, yet your data ends at 5. Either you normalize it from 0-1, or add a new parameter to scale your data (that's what I did).
2. Your x values are too large and you're getting overflows. I divided them by 1E5 and I got an appropriate fit.

Here's my code. Note I changed the import of matplotlib. Other than that, it's mostly the same.

import numpy as np
import matplotlib.pyplot as plt # Changed here
from scipy.optimize import curve_fit

def sigmoid(x, a, b, c): # Added param c
    return c / (1.0 + np.exp(-a*(x-b)))

# Made xdata into an array
xdata = np.array([ 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]

xdata = xdata / 1E5  # scaled it down
         
popt, pcov = curve_fit(sigmoid, xdata, ydata, method='dogbox', bounds=([min(ydata), min(xdata), 0],[max(ydata), max(xdata), 10])) # added some bounds for the `c` parameter
print(popt)

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

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

Here's the result: popt=[4.40402529 0.93238322 5.13341602],