How can I reduce the artifacts generated during "Thin Plate Spline" interpolation in Python?

Question

At the Top "right", there is the 2D-density plot of the recorded data (actual), fewer in number. On the top-left is the interpolated data (thin-plate), i.e. larger in number. Compared to the Right (actual) data-set we can observe that "Interpolated" data is showing some "Butterfly" kind of artifact which is not desired and expected. I have tried both my own written-code and in-built python module to perform "Thin-Plate-spline" interpolation and both the results are showing same "Butterfly" artifact.

Moreover, I am looking for "smooth" 3D-surface as demonstrated in 2 figures at the bottom.

Code 1 => Density-Plot

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

###########################################################
#np.random.seed(1977)
#x, y, z = np.random.random((3, 10))
###########################################################
data = np.loadtxt('SIM24_List.dat',unpack=True)
data1 = np.loadtxt('data99.dat',unpack=True)

x = (data[0,:])/100 
y = (data[1,:])/100
E1 = data1[0,:]
###########################################################

interp = scipy.interpolate.Rbf(x, y, E1, function='thin_plate')

#yi, xi = np.mgrid[0:1:100j, 0:1:100j]

###########################################################
X = np.linspace(-300, 300, 1200)
Y = np.linspace(-300, 300, 1200)
xi, yi = np.meshgrid(X, Y)
###########################################################

Ei = interp(xi, yi)

###########################################################
plt.plot(x, y, 'ko')
plt.imshow(Ei, extent=[0, 1, 1, 0], cmap='gist_earth')
plt.colorbar()
###########################################################
plt.show()

Code 2 => 3D-Surface Plot

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

#############################################################
from mpl_toolkits.mplot3d import axes3d  #Axes3d
import matplotlib.cm as cm
from pylab import rcParams
rcParams ['figure.figsize'] = 5, 5
#########################################################################################################

fig = plt.figure()
axes = fig.gca(projection='3d') #gca = get current axis

###########################################################
#np.random.seed(1977)
#x, y, z = np.random.random((3, 10))
###########################################################
data = np.loadtxt('SIM24_List.dat',unpack=True)
data1 = np.loadtxt('data99.dat',unpack=True)

x = (data[0,:])/100 
y = (data[1,:])/100
E1 = data1[0,:]
###########################################################

interp = scipy.interpolate.Rbf(x, y, E1, function='thin_plate')

#yi, xi = np.mgrid[0:1:100j, 0:1:100j]

###########################################################
X = np.linspace(-300, 300, 1200)
Y = np.linspace(-300, 300, 1200)
xi, yi = np.meshgrid(X, Y)
###########################################################

Ei = interp(xi, yi)

###########################################################
#p = axes.plot_surface(xi,yi,Ei, rstride=4, cstride=4, cmap=cm.RdBu, linewidth=0, antialiased=False)
p = axes.plot_surface(xi, yi, Ei,cmap='viridis', edgecolor='none')
fig.colorbar(p, shrink=0.5)
plt.tight_layout();
fig.savefig("surface.pdf")
###########################################################
plt.show()

Code 3 => My own code - 3D

import numpy as np
import math as mt


import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d  #Axes3d
import matplotlib.cm as cm
from pylab import rcParams
rcParams ['figure.figsize'] = 5, 5

#########################################################################################################

fig = plt.figure()
axes = fig.gca(projection='3d') #gca = get current axis


#########################################################################################################
print(" ")
################################################################

#data = np.loadtxt('data1.dat',unpack=True)
#x = data[0,:] 
#y = data[1,:]
#E1 = data[2,:]

################################################################

data = np.loadtxt('SIM24_List.dat',unpack=True)
data1 = np.loadtxt('data99.dat',unpack=True)

x = (data[0,16:])/100 
y = (data[1,16:])/100
E1 = data1[0,16:]

################################################################
################################################################
n1 = len(x)
#print(" n1 : ", n1)
r = np.stack((x, y), axis=-1)
print(r)


N = n1 + 3

o = np.zeros(3)


#################################################################

#################################################################

V1 = np.hstack((E1, o))
V = V1.reshape(N, 1)

##########################################################################################################


C = np.ones(n1)

P = np.stack((C, x, y), axis=-1)
P_T = P.transpose()  # Transpose of P

#print(" P : ")
#print(P)

#print(" P_T : ")
#print(P_T)


#########################################################################################################

O1 = np.zeros(9)
O = O1.reshape((3,3))

#print(" O : ")
#print(O)

#########################################################################################################

K = np.zeros((n1,n1))

for i in range(n1):
    for j in range(n1):

        if(i == j):
            K[i][j] = 0.0
        else:
            R = r[i] - r[j]
            s = np.sqrt( (R[1]*R[1]) + (R[0]*R[0]) )
#           print(s)
            K[i][j] =  (s*s) * (np.log(s)) 

#print(" U[i][j]")
#print(K)

##########################################################################################################

#L1 = np.block([K, P])
L1 = np.hstack((K, P))
#L2 = np.block([P_T, O])
L2 = np.hstack((P_T, O))
#L = np.block([[L1],[L2]])
L = np.vstack((L1,L2))
#print(" L : ")
#print(L)

LI = np.linalg.inv(L)

###########################################################################################################

W = np.dot(LI, V)   # Weights

print(" ")
#print(" W : ")
#print(W)

###########################################################################################################

def f(x, y, r, N, W):

    f = W[N-3] + (W[N-2])*x + (W[N-1])*y
 
    for i in range(N-3):

        r1 = np.sqrt(((r[i][0] - x)**2) + ((r[i][1] - y)**2))
#       r1 = (((r[i][0] - x)**2) + ((r[i][1] - y)**2))
        if(r1==0.0):
            U1 = 50     #IMPORATANT
        else:
            U1 = (r1*r1) * (np.log(r1))

                
        f = f + W[i]*U1
    
    if(f < 0.0):
        f = 0.0

    return f


#f9 = f(35, 35, r, N, W)  # value of f9 is in good approximation with the value achieved from simulation
#print(f9)

## up to this point CODE is working well


print("_______________________________________________________________________________________")
print(" ")


#F = f(xv, yv, x, y, N, W) # This is not working 



x1 = np.linspace(-200, 200, 1000)
y1 = np.linspace(-200, 200, 1000)
xv, yv = np.meshgrid(x1,y1)

f_xy = np.zeros((1000, 1000))

for i in range(1000):
    for j in range(1000):
        f_xy[i][j] = f(x1[i], y1[j], r, N, W)


#############################################################################################################################

p = axes.plot_surface(xv,yv,f_xy, rstride=4, cstride=4, cmap=cm.RdBu, linewidth=0, antialiased=False)
##p = axes.plot_surface(xv, yv, f_xy,cmap='viridis', edgecolor='none')
fig.colorbar(p, shrink=0.5)
plt.tight_layout();
fig.savefig("surface.pdf")
plt.show()