Biexponential fit doesn't match with data

I have a set of data for which I am trying to fit a biexponential function. So far, I have done

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import csv
from pylab import genfromtxt
import scipy

twofifty = genfromtxt("t vs delta E (250C).txt");

t = twofifty[:, 0]
de = twofifty[:, 1]*10e3

def func(t, f_1, tau_1, f_2, tau_2):
    return f_1*np.exp(-t/tau_1)+f_2*np.exp(-t/tau_2)

guess = (0.7, 100, 0.3, 200)

popt, pcov = curve_fit(func, t, de, guess)
print(popt)
plt.plot(t, de, label="data")
plt.plot(t, func(de, *popt), label="fit")
plt.legend()
plt.show()

But the graph comes out as shown below.

I'm not sure what I need to change to get the correct results. I have changed the scale of the y axis since it was in the order of 10^4 and included some guess values which allowed me to get the following results, [-1.46219808e+02 -4.99772483e+05 2.42937404e+01 -1.94190603e+01], but I know they are not the correct values as the fit is not right and they are negative values. As you can tell, from the biexponential equation I am trying to extract the amplitudes (f_1 and f_2) and decay time (tau_1 and tau_2). I know for the amplitudes the sum has to be 1, hence my predictions of 0.7 and 0.3, but as I don't really know what tau will be, only that they will be of the order of 10^-1 or 10^-2. How can I smooth out the fit line, because as it's evident, it doesn't actually fit the data? What can I change to make it fit the data?

When doing the change of plt.plot(t, func(de, *popt), label="fit") to plt.plot(t, func(t, *popt), label="fit"), as recommended, the plot comes out to be

This smooths out the curve, but still doesn't seem to match/fit the data.

Your are facing a simple problem flat baseline when t <= 0 cannot be captured by your model which will have to grow asymptotically.

Censuring data when t <= 0, we got:

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize, stats, special

x, y = np.genfromtxt("https://pastebin.com/raw/h2HnULEz").T
q = x >= 0
x = x[q]
y = y[q]

def model(t, f1, tau1, f2, tau2):
    return f1 * np.exp(- t / tau1) + f2 * np.exp(- t / tau2)

popt, pcov = optimize.curve_fit(model, x, y, maxfev=10000)

By removing the flat baseline before and forcing more function evaluations we get a bit more fitness:

But this is clearly not satisfying and certainly does not capture the initial flat baseline.

On the other hand your data strongly looks like a negative version of a Log Normal distribution. Writing the model and fitting:

def model2(t, A, sigma, loc, scale):
    law = stats.lognorm(sigma, loc=loc, scale=scale)
    return - A * law.pdf(t)

popt2, pcov2 = optimize.curve_fit(model2, x, y, p0=[0.1, 0.5, 1, 10])

Gives in comparison a better fit:

If we keep the whole data, without any precaution but an educated guess we got:

popt, pcov = optimize.curve_fit(model, x, y, p0=[30, 5, -30, 5])

If we clamp the asymptotic behaviour around t < 2. we get almost a good fit, but not better than the Log Normal:

def model(t, f1, tau1, f2, tau2):
    y = f1 * np.exp(- t / tau1) + f2 * np.exp(- t / tau2)
    y[t < 2] = 0.
    return y

This is strong evidence that the main problem with your fit is the asymptotic behaviour of the model that cannot be captured without clamping your model or censoring your dataset.