How to plot a curve fit over a given range in matplotlib

I'm trying to replicate the following plot in python, specifically the dotted "Partial Superparticle Ionisation" curve fit. Here's what my current fit looks like. I can't seem to figure out how to get it to start at 0 and then level off before 1.

My Python code is as follows:

import numpy as np
import numpy as np
import scipy.optimize as opt;
import matplotlib.pyplot as plt 

x = np.array([0.3, .6, .7, .8, 1.3, 1.5, 3, 4, 5])
y = np.array([0, .1, .3, .55, .75, .9, 1, 1, 1])

def func(x, a, b, c):
     return a * np.exp(-b * x) + c

# The actual curve fitting happens here
optimizedParameters, pcov = opt.curve_fit(func, x, y);

# Use the optimized parameters to plot the best fit
plt.plot(x, func(x, *optimizedParameters), label="fit");

# Plotting the Graph
plt.ylim(0, 1.1)
plt.xlim(0, 5)
plt.step(x, y)
plt.xlabel("Time (ns)")
plt.ylabel("$n_H^+ / n_H$")
plt.show()

To allow this to start at 0,0 you first need to change the func to:

def func(x, a, b):
    return a * (1 - np.exp(-b * x))

Then you need to add the bounds parameter to the opt.curve_fit function:

optimizedParameters, pcov = opt.curve_fit(func, x, y, bounds=([0, 0], [1, 5]))

Finally you need to generate x values for the fitted curve:

# Generate x values for the fitted curve
x_fit = np.linspace(0, 5, 500)
y_fit = func(x_fit, *optimizedParameters)

The full code script is as follows:

import numpy as np
import numpy as np
import scipy.optimize as opt;
import matplotlib.pyplot as plt 

x = np.array([0, 0.3, .7, .8, 1.3, 1.5, 3, 5])
y = np.array([0, .1, .3, .55, .75, .9, 1, 1])

def func(x, a, b):
    return a * (1 - np.exp(-b * x))

optimizedParameters, pcov = opt.curve_fit(func, x, y, bounds=([0, 0], [1, 5]))

# Generate x values for the fitted curve
x_fit = np.linspace(0, 5, 500)
y_fit = func(x_fit, *optimizedParameters)

# Plotting the Graph
plt.figure()
plt.step(x, y, where='post', label="Whole Superparticle Ionisation", color='blue')
plt.plot(x_fit, y_fit, 'k--', label="Partial Superparticle Ionisation")
plt.ylim(0, 1.1)
plt.xlim(0, 5)
plt.xlabel("Time (ns)")
plt.ylabel("$n_{H^+} / n_H$")
plt.legend()
plt.show()