Some code which computes Fourier coefficients for the periodic Fourier series. Then it plots both the function and the partial sum Fourier approximation.

import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

# -- If you want to give the depth when you execute, then use this code:
# from sys import argv
# script, depth_provided, filename = argv
# depth = int(depth_provided)
# -- otherwise, set the depth and filename here:
depth = 5
filename = "Example.pdf"

# -- the function being approximated --
u = lambda x: x**3+1

# -- the domain is [-L,L] --
L = np.pi


a = range(0,depth+1,1)
b = range(0,depth+1,1)
	
# -- define the normalized sine and cosine functions -- 
# -- (Note: these can be replaced by any orthonormal functions!) --	

def normal_sine(k): return lambda x: np.sqrt(1./L)*np.sin(k*x)

def normal_cosine(k): return lambda x: np.sqrt(1./L)*np.cos(k*x)

normal_const = lambda x: 1./np.sqrt(2.*L) + 0.*x

# -- function operations --
def mult(u,v): return lambda x: u(x)*v(x)
def add(u,v): return lambda x: u(x) + v(x)

# -- compute a[0], the constant term in the series --
a[0] = quad( mult(u,normal_const) ,-L,L)[0]

# -- compute other coefficients --
for k in range(1,depth+1,1):
	a[k] = quad( mult(u,normal_cosine(k)) ,-L,L)[0]

for k in range(1,depth+1,1):
	b[k] = quad(mult(u,normal_sine(k)),-L,L)[0]
	
# -- define the approximate function --	
def fourier_approx(x):
	value = a[0]*normal_const(x)
	for k in range(1,depth+1,1):
		value = value + a[k]*normal_cosine(k)(x) + b[k]*normal_sine(k)(x)
	return(value)	

# -- plot the original function and the Fourier approximate --
x_values = np.arange(-L,L,0.001)

plt.plot(x_values, fourier_approx(x_values),label="FS (depth %r)" %depth)
plt.plot(x_values, u(x_values),label=r"$u(x)$")

plt.legend(loc="upper left")

plt.savefig(filename)
# -- if you want to actually see the plot in "real time" --
#plt.show()