I am trying to perform 2D bicubic spline interpolation following this article (pg 15 sec 6.2 Bicubic Interpolation)

My test -- I have a gaussian function defined on the interval $ x,y \in [-1,1] $ and I want to use the bicubic interpolation method to get a reasonable interpolation of the function.

This is my implementation so far -

import numpy as np  
from numdifftools import Hessian, Gradient

def gaussian(pt):
    x, y = pt
    return np.exp(-x**2 - y**2)

A_inv = np.array([
    [1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],
    [-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0],
    [2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0],
    [0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0],
    [0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0],
    [0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0],
    [-3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0],
    [0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0],
    [9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1],
    [-6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1],
    [2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0],
    [0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0],
    [-6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1],
    [4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1]
])


def get_coefficient_tensor(func_array):
    alpha_array = A_inv @ func_array
    return np.array([
        [alpha_array[0], alpha_array[1], alpha_array[2], alpha_array[3]],
        [alpha_array[4], alpha_array[5], alpha_array[6], alpha_array[7]],
        [alpha_array[8], alpha_array[9], alpha_array[10], alpha_array[11]],
        [alpha_array[12], alpha_array[13], alpha_array[14], alpha_array[15]]
    ]).T

def constructed_function(x, y, x_range, y_range, func_array):
    ''' 
    Returns the value of the constructed interpolation function at the point (x,y)

    Parameters:
    x : float
        x coordinate of the point
    y : float
        y coordinate of the point
    x_range : list / touple of 2 floats
        range of the x coordinate of the current square grid being considered
    y_range : list / touple of 2 floats
        range of the y coordinate of the current square grid being considered
    func_array : numpy array of size 16
        array containing the user defined quantities - 16x1 dimensional array
            - first 4 elements are the values of the function at the 4 corners of the square
            - next 8 elements are the values of the partial derivatives in the two directions at the 4 corners of the square (susceptibilities)
            - last 4 elements are the values of the mixed partial derivatives at the 4 corners of the square (I don't know how to interpret this yet)
    '''
    t = (x - x_range[0]) / (x_range[1] - x_range[0])
    u = (y - y_range[0]) / (y_range[1] - y_range[0])
    coefficient_tensor = get_coefficient_tensor(func_array)
    return np.sum([
        coefficient_tensor[i,j] * t**i * u**j for i in range(4) for j in range(4)
    ])

def get_func_array(func, x_range, y_range):
    # function values at the 4 corners
    f_00 = func([x_range[0], y_range[0]])   
    f_01 = func([x_range[0], y_range[1]])
    f_10 = func([x_range[1], y_range[0]])
    f_11 = func([x_range[1], y_range[1]])

    # partial derivatives at the 4 corners
    f_x_00, f_y_00 = Gradient(func)([x_range[0], y_range[0]])
    f_x_01, f_y_01 = Gradient(func)([x_range[0], y_range[1]])
    f_x_10, f_y_10 = Gradient(func)([x_range[1], y_range[0]])
    f_x_11, f_y_11 = Gradient(func)([x_range[1], y_range[1]])

    # cross partial derivatives at the 4 corners
    f_xy_00 = Hessian(func)([x_range[0],y_range[0]])[0,1] 
    f_xy_01 = Hessian(func)([x_range[0],y_range[1]])[0,1] 
    f_xy_10 = Hessian(func)([x_range[1],y_range[0]])[0,1] 
    f_xy_11 = Hessian(func)([x_range[1],y_range[1]])[0,1] 
    return np.array([f_00,f_01,f_10,f_11,f_x_00,f_x_01,f_x_10,f_x_11,f_y_00,f_y_01,f_y_10,f_y_11,f_xy_00,f_xy_01,f_xy_10,f_xy_11]).T
                                      

x_range = [-1,1]
y_range = [-1,1]

xs = np.linspace(x_range[0], x_range[1], 100)
ys = np.linspace(y_range[0], y_range[1], 100)

X, Y = np.meshgrid(xs, ys)

func_array = get_func_array(gaussian, x_range, y_range)

Z = np.array([
    constructed_function(x, y, x_range, y_range, func_array) for x in xs for y in ys
]).reshape(X.shape)

Z_actual = np.array([
    gaussian([x,y]) for x in xs for y in ys
]).reshape(X.shape)

When I visualize the Z_interpolated (green), it is very different from Z_actual (red) -- not only are the derivatives the wrong magnitude at the corners, they are also the wrong sign.

I want the interpolation scheme to better capture the derivatives at the four points I have specified it.

Is this a limitation of the bicubic spline or is there an error in the implementation?

Should I include more grid points? I feel like since it is a piecewise interpolation scheme should it really improve things?