# %% Import of packages:

import matplotlib.pyplot as plt
import numpy as np

# %% Functions calculating alternative num derivatives:

def fun_f(x):
    f = x**3 + 1    
    return f

# %% Generating x and y:

x_start = 0
x_stop = 1
dx = 0.1
N = int((x_stop - x_start)/dx + 1)  # Number of samples
x_array = np.linspace(x_start, x_stop, N)
# Not included if y is an existing data series:
y_array = fun_f(x_array)

# %% Preallocation of arrays:
# np.nan (not-a-number) makes plotting code simpler.

dy_dx_center_array = np.zeros(N) + np.nan

# %% Calculation of the derivative in a for loop:

for k in range(1, N-1):
    dy_dx_center_array[k] = (y_array[k+1] - y_array[k-1])/(2*dx)
        
# %% Continuous dy_dx for comparison with center dy_dx:

dy_dx_cont_array = 3*x_array**2

x_cont_start = 0
x_cont_stop = 1
dx_cont = 0.001
N_cont = int((x_stop - x_start)/dx_cont + 1)
x_cont_array = np.linspace(x_start, x_stop, N_cont)

y_cont_array = x_cont_array**3 + 1
dy_dx_cont_array = 3*x_cont_array**2

# %% Plotting:

plt.close('all')  # Closes all figures before plotting
plt.figure(num=1, figsize=(24/2.54, 18/2.54)) # Inch

plt.subplot(2, 1, 1)
plt.plot(x_array, y_array, 'or', 
         label='y')
plt.plot(x_cont_array, y_cont_array, 'b-', 
         label='y_cont')
plt.legend()
plt.grid(which='both', color='grey')


plt.subplot(2, 1, 2)
plt.plot(x_array, dy_dx_center_array, 'or', 
         label='dy_dx_center')
plt.plot(x_cont_array, dy_dx_cont_array, 'b-', 
         label='dy_dx_cont')
plt.grid(which='both', color='grey')
plt.legend()
plt.xlabel('x')

# plt.savefig('num_diff_center.pdf')  # pdf
plt.show()