diff --git a/examples/example.py b/examples/example_dirichlet.py
similarity index 100%
rename from examples/example.py
rename to examples/example_dirichlet.py
diff --git a/examples/example_neumann.py b/examples/example_neumann.py
new file mode 100644
index 0000000000000000000000000000000000000000..826dba14b6b9c5c70e436bec3e142c7b431aa672
--- /dev/null
+++ b/examples/example_neumann.py
@@ -0,0 +1,92 @@
+# Example file to demonstrate the use of the VKOGA-PDE code for
+# a PDE that uses both Dirichlet and Neumann boundary conditions.
+# We use \Omega = [0, 2]^2, L = -LaPlace and a Wendland kernel.
+
+
+
+
+# Some imports
+import numpy as np
+from vkoga_pde.kernels_PDE import wendland33_laplace
+from vkoga_pde.vkoga_PDE import VKOGA_PDE
+from matplotlib import pyplot as plt
+import math
+
+np.random.seed(1)
+
+## Set number of points that will be sampled
+N1 = int(5e3)
+N2 = int(1e2)
+N3 = int(1e2)
+
+## Create some data to describe domain (X1), Dirichlet (X2) and Neumann (X3) boundary
+dim = 2
+
+X1 = np.random.rand(N1, dim)
+
+X2 = np.random.rand(2*N2, dim)
+X2[0:N2, 1] = 0
+X2[N2:, 1] = 1
+
+X3 = np.random.rand(2*N3, dim)
+X3[:N3, 0] = 1
+X3[N3:, 0] = 0
+
+n_X3 = np.zeros((2*N3, dim))
+n_X3[:N3, 0] = 1
+n_X3[:N3, 0] = 0
+
+n_X3[N3:, 0] = -1
+n_X3[N3:, 1] = 0
+
+# Rescale the domain a bit: Kernel width parameters are not implemented correctly, thus just scale the domain
+scale_factor = 2
+X1 = scale_factor * X1
+X2 = scale_factor * X2
+X3 = scale_factor * X3
+
+# Visualize the domain
+plt.figure(1)
+plt.clf()
+plt.plot(X1[:, 0], X1[:, 1], 'm.')
+plt.plot(X2[:, 0], X2[:, 1], 'r.')
+plt.plot(X3[:, 0], X3[:, 1], 'b.')
+
+
+## Define functions u (solution of PDE), f (-Laplace u = f) and g (for Neumann)
+u = lambda x: 20 * x[:, [1]] * x[:, [0]] + 1
+f = lambda x: np.zeros((x.shape[0], 1))
+g = lambda x: 20 * x[:, [1]] * np.sign(x[:, [0]] - .5)
+
+# Define right hand side values
+y1 = f(X1)      # Application of -Laplace
+y2 = u(X2)      # Dirichlet boundary, given by solution
+y3 = g(X3)      # Neumann boundary
+
+
+## Initialize and run models
+kernel = wendland33_laplace(dim=2)        # or also try this une
+maxIter = 100
+
+model_pde = VKOGA_PDE(kernel=kernel, greedy_type='f_greedy')
+_ = model_pde.fit(X1, y1, X2, y2, X3=X3, y3=y3, n_X3=n_X3, maxIter=maxIter)#, y1_sol=u(X1))
+
+
+
+## Plot the solution
+N_plots = 1000              # don't plot too many points, otherwise plot is slow
+
+fig = plt.figure(105)
+ax = fig.add_subplot(projection='3d')
+ax.scatter(X1[:N_plots, 0], X1[:N_plots, 1], model_pde.predict_s(X1[:N_plots, :]))
+ax.set_xlabel('X Label')
+ax.set_ylabel('Y Label')
+ax.set_zlabel('Z Label')
+#ax.set_zlim(-2,2)
+plt.show()
+
+
+
+
+
+