Why is symfit giving me terrible results when I try to fit fourier series to my data?

Question

I am trying to find a fourier series that fits this plot:

First I used this website and I found the following fourier coefficients

a = [0.983, 0.292, 0.110, 0.006, -0.022, -0.018, 0.006]
b = [0.000, 0.246, 0.002, -0.021, 0.017, 0.046, 0.019]
omega = 0.045 

which produced the following result

Which is pretty good, but then I decided to try to code my own curve fitter in python, and found that the library symfit seemed to be the best option. I even repurposed one example found in the documentation while changing very little.

And sure enough when I run the code the fit seems to be pretty good in the graphs made with the fit object directly:

plt.plot(xdata, fit.model(x=xdata, **fit_result.params).y, color='green', ls=':')

But then I decided to use the fourier coefficients directly and not rely on the fit object to try to reproduce the function. This is how the code prints those coefficients:

Parameter Value        Standard Deviation
a0        -1.206104e+03 nan
a1        1.799223e+03 1.087713e+05
a10       2.393892e+00 8.023119e+02
a2        -8.465307e+02 2.133670e+05
a3        6.798590e+02 nan
a4        -9.287399e+02 nan
a5        7.388156e+02 nan
a6        -2.099387e+02 1.376980e+05
a7        -9.838123e+01 9.799211e+04
a8        9.821148e+01 3.739932e+04
a9        -2.827002e+01 8.152094e+03
b1        -6.084419e+02 2.480938e+05
b10       4.586827e+01 3.396681e+02
b2        9.311030e+02 2.178771e+05
b3        -8.877734e+02 2.406878e+05
b4        8.810145e+02 4.462669e+05
b5        -1.288586e+03 4.376938e+05
b6        1.734395e+03 2.730319e+05
b7        -1.583775e+03 1.137506e+05
b8        9.115344e+02 3.084382e+04
b9        -3.045882e+02 4.905459e+03
w         -1.370387e-02 nan

The standard deviations were a huge red flag, they are absurdly huge, and sure enough, this is the result

Which seems like a bad joke.

That was done using 10 coefficients, if I only use six, I get something barely acceptable.

The worst part is that the graphs made directly with the fit object must mean that the code did produce a good fit, but for some reason the fourier coefficients it gives me are just lies.

Here is my code.

First, this is how I fit the functions

from symfit import parameters, variables, sin, cos, Fit
import numpy as np
import matplotlib.pyplot as plt
import math as m

def fourier_series(x, f, n=0):
    a0,*cos_a = parameters(','.join(['a{}'.format(i) for i in range(0, n + 1)]))
    sin_b = parameters(','.join(['b{}'.format(i) for i in range(1, n + 1)]))
    series = a0 + sum(ai * cos(i * f * x) + bi * sin(i * f * x)
                     for i, (ai, bi) in enumerate(zip(cos_a, sin_b), start=1))
    return series

x, y = variables('x, y')
w, = parameters('w')
model_dict = {y: fourier_series(x, f=w, n=10)}
print(model_dict)

#This is the graph I'm trying to fit turned into data points
xlist=[381.854,385.274,387.55,388.973,390.656,392.671,394.281,396.37,398.014,402.047,408.288,411.885,414.452,418.099,\
       423.493,426.901,428.014,430.068,433.114,438.699,445.023,448.562,453.307,458.014,463.663,467.877,474.536,\
       479.384,485.409,489.247,496.283,501.164,508.191,514.315,520.1]
ylist=[0.02816,0.09827,0.18838,0.2948,0.36316,0.46243,0.53795,0.60116,0.64162,0.68361,0.71676,0.79285,0.87861,0.94579,\
       0.99422,0.95307,0.89017,0.85549,0.79285,0.76301,0.7273,0.65896,0.59621,0.53757,0.47969,0.43353,0.38501,\
       0.33526,0.28305,0.23121,0.18838,0.14451,0.11555,0.08092,0.07185]

xdata = np.array(xlist)
ydata = np.array(ylist)

# Define a Fit object for this model and data
fit = Fit(model_dict, x=xdata, y=ydata)
fit_result = fit.execute()
print(fit_result)

# Plot the result
plt.plot(xdata, ydata)
plt.plot(xdata, fit.model(x=xdata, **fit_result.params).y, color='green', ls=':')

This is how I take those coefficients and try to reproduce the graph. I have different sets of a b and omega data to try to reproduce other fourier series. The one for the step wedge (the example in the documentation) is the only one that works.

from symfit import parameters, variables, sin, cos, Fit
import numpy as np
import matplotlib.pyplot as plt
import math as m
from random import random

def f(x,a,b,omega):
    output = a[0]/2
    for k in range(1,len(a)-1):
        output = output + a[k]*m.cos(k*omega*x) + b[k]*m.sin(k*omega*x)
    return output

def monteCarlo(min,max,a,b,omega):
    x = min + random()*(max-min)
    return (x,f(x,a,b,omega))

xdata=[]
ydata=[]

#for step function between -3 and 3
#a=[5.000000e-01,1.106075e-10,2.364473e-11,1.697959e-10,2.822410e-11,3.860925e-11,-1.522386e-10]
#b=[0,6.267589e-01,2.020945e-02,1.846406e-01,3.623074e-02,8.726419e-02,4.518721e-02]
#omega=8.615115e-01

#for our spectrum with n=6
a = [-1.755069e+02,-1.097051e+02,1.571135e+02,6.491389e+01,1.220463e+02,1.778165e+02,4.564577e+01]
b = [0,-5.523502e+01,1.517397e+02,1.325198e+02,-1.015198e+02,-3.255781e+01,2.495038e+01]
omega = -5.353751e-03

#for our spectrum with n=10
#a = [-1.206104e+03,1.799223e+03,-8.465307e+02,6.798590e+02,-9.287399e+02,7.388156e+02,-2.099387e+02,-9.838123e+01,\
#    9.821148e+01,-2.827002e+01,2.393892e+00]
#b = [0,-6.084419e+02,9.311030e+02,-8.877734e+02,8.810145e+02,-1.288586e+03,1.734395e+03,-1.583775e+03,9.115344e+02,\
#     -3.045882e+02,4.586827e+01]
#omega = -1.370387e-02

#a = [0.983, 0.292, 0.110, 0.006, -0.022, -0.018, 0.006]
#b = [0.000, 0.246, 0.002, -0.021, 0.017, 0.046, 0.019]
#omega = 0.045

for n in range(0,1000):
    coor = monteCarlo(375,525,a,b,omega)
    xdata.append(coor[0])
    ydata.append(coor[1])

plt.plot(xdata, ydata, 'bo')