diff --git a/pybie2d/boundaries/global_smooth_boundary/global_smooth_boundary.py b/pybie2d/boundaries/global_smooth_boundary/global_smooth_boundary.py
index b2e594e..3eefc1d 100644
--- a/pybie2d/boundaries/global_smooth_boundary/global_smooth_boundary.py
+++ b/pybie2d/boundaries/global_smooth_boundary/global_smooth_boundary.py
@@ -1,5 +1,6 @@
 import numpy as np
 import scipy as sp
+import matplotlib as mpl
 import scipy.signal
 import warnings
 import os
@@ -74,6 +75,7 @@ def __init__(self, x=None, y=None, c=None):
         self.max_h = np.max(self.weights)
         self.area = self.dt*np.sum(self.x*self.cp.imag)
         self.perimeter = self.dt*np.sum(self.speed)
+        self.poly= mpl.path.Path(np.column_stack([self.x, self.y]))
 
         self.defined_modules =  [   'Laplace_SLP_Self_Kress',
                                     'Stokes_SLP_Self_Kress',
@@ -193,15 +195,49 @@ def Get_Close_Corrector(self, kernel, *args, **kwargs):
     #########################
     #### Private Methods ####
     #########################
-
-    def _test_inside_point(self, candidate, eps=1e-10):
+    
+    # THIS ROUTINE IS NOT RELIABLE
+    #def _test_inside_point(self, candidate, eps=1e-10):
+    #    """
+    #    Test whether the provided or generated inside point is acceptable
+    #    returns True if the point is okay, False if its not
+    #    """
+    #    test_value = np.sum(self.complex_weights/(self.c-candidate))
+    #    return np.abs(test_value - 2.0j*np.pi) < eps
+    ## end _test_inside_point function
+ 
+    # THIS ROUTINE IS RELIABLE AND INDEPENDENT OF MATPLOTLIB
+    #def _test_inside_point(self, candidate, eps=1e-12):
+    #    """
+    #    Test whether the provided or generated inside point is acceptable
+    #    returns True if the point is okay, False if its not
+    #    Algorithm: point in polygon (the points in self.c are interpreted as closed polygon),
+    #    the sum of angles between the vectors from the point (candidate) to two consecutive points on the boundary #
+    #     polygon. the sum is equal 2pi only if the point lies inside. Safe way to compute the angle is using arctan2 of #(a x b)/(a.b).
+    #     """
+    #     a=self.c-candidate
+    #     if (np.any(np.abs(a)<eps)):
+    #       return False # closer than eps to a boundary point, so labelled not inside
+    #     # angle between two "vectors" expressed as a complex number: a[i]=c[i]-cpoint,b=a[i+1]
+    #     b=np.roll(a,-1,axis=0)
+    #     axb=(a.real*b.imag-a.imag*b.real)
+    #     adotb=a.real*b.real+a.imag*b.imag
+    #     # compute sum of angles
+    #     test_value=np.sum(np.arctan2(axb,adotb)) 
+    #     return np.abs(test_value - 2.0*np.pi) < eps
+    ##  end _test_inside_point function
+ 
+    # THIS ROUTINE IS RELIABLE, DEPENDS OF MATPLOTLIB (same use as in `find_interior_points`) 
+    def _test_inside_point(self, candidate):
         """
         Test whether the provided or generated inside point is acceptable
-        returns True if the point is okay, False if its not
+        returns True if the point (=candidate) is okay, False if its not
+        Algorithm: 
+        approximate the boundary points as a closed polygon, which is
+        initialized in self.poly using matplotlib.path. Then use the contains_point method
+        to check if the point lies inside the polygon (analogous to `find_interior_points`)
         """
-        test_value = np.sum(self.complex_weights/(self.c-candidate))
-        return np.abs(test_value - 2.0j*np.pi) < eps
-    # end _test_inside_point function
+        return self.poly.contains_point([candidate.real,candidate.imag])
 
 def arr_check(x):
     if type(x) != np.ndarray:
diff --git a/tests/test_global_smooth_boundary.ipynb b/tests/test_global_smooth_boundary.ipynb
new file mode 100644
index 0000000..993f347
--- /dev/null
+++ b/tests/test_global_smooth_boundary.ipynb
@@ -0,0 +1,105 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%reload_ext autoreload\n",
+    "%autoreload 2\n",
+    "import numpy as np\n",
+    "import pybie2d\n",
+    "import matplotlib as mpl\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "star = pybie2d.misc.curve_descriptions.star\n",
+    "GSB = pybie2d.boundaries.global_smooth_boundary.global_smooth_boundary.Global_Smooth_Boundary\n",
+    "\n",
+    "def check_GSB_test_inside_points(NBfactor):\n",
+    "    '''\n",
+    "    Test the boundary._test_inside_point method, using  points that are known to be inside/outside a star-shaped closed globally-smooth boundary\n",
+    "    '''\n",
+    "    \n",
+    "    NB0=3*5\n",
+    "    NB=NB0*NBfactor\n",
+    "    boundary = GSB(c=((star(NB,x=1,y=1,r=1,a=0.6,f=7))))\n",
+    "    \n",
+    "    print(f\"NB={NB}\")\n",
+    "    #two points on the boundary\n",
+    "    pc = boundary.c[[1*NB//5,2*NB//3]]\n",
+    "    eps=1e-10\n",
+    "    t=np.linspace(eps,1-eps,30)\n",
+    "    inner_pts=pc[0]+(pc[1]-pc[0])*t-0.5*np.sin(np.pi*t)*(pc[1].real-pc[0].real)\n",
+    "    outer_pts=np.concatenate((pc[1]+(pc[1]-pc[0])*t+0.5*np.sin(np.pi*t)*(pc[1].real-pc[0].real),\n",
+    "                              pc[0]-(pc[1]-pc[0])*t+0.5*np.sin(np.pi*t)*(pc[1].real-pc[0].real)))\n",
+    "\n",
+    "\n",
+    "    \n",
+    "    check_inpts=np.array([False]*len(inner_pts))\n",
+    "    for i,inpt in enumerate(inner_pts):\n",
+    "        check_inpts[i]=boundary._test_inside_point(inpt)\n",
+    "    \n",
+    "    if not np.all(check_inpts): \n",
+    "        print(f'test_inside_point failed for {np.sum((~check_inpts)*1)} /{30} points that are inside the boundary')\n",
+    "    failed_inpts=inner_pts[~check_inpts]\n",
+    "    success_inpts=inner_pts[check_inpts]\n",
+    "\n",
+    "   \n",
+    "    check_outpts=np.array([False]*len(outer_pts))    \n",
+    "    \n",
+    "    for i,outpt in enumerate(outer_pts):\n",
+    "        check_outpts[i]=~boundary._test_inside_point(outpt)\n",
+    "    \n",
+    "    if not np.all(check_outpts): print(f'test_inside_point failed for {np.sum((~check_outpts)*1)} /{60} points that are outside the boundary')\n",
+    "\n",
+    "    failed_outpts=outer_pts[~check_outpts]\n",
+    "    success_outpts=outer_pts[check_outpts]\n",
+    "    failed_pts=np.concatenate((failed_inpts,failed_outpts))\n",
+    "    success_pts=np.concatenate((success_inpts,success_outpts))\n",
+    "    return boundary.c,failed_pts,success_pts\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "pts,failed_pts,success_pts=check_GSB_test_inside_points(10)\n",
+    "fig, ax = plt.subplots(1,1,figsize=(12,12))\n",
+    "ax.set_aspect('equal')\n",
+    "\n",
+    "ax.plot(pts.real,pts.imag,'k-')\n",
+    "ax.plot(failed_pts.real,failed_pts.imag,'rx')\n",
+    "ax.plot(success_pts.real,success_pts.imag,'bx') \n",
+    "\n",
+    "assert(len(failed_pts)==0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}