diff --git a/INSTALL.MacOSX b/INSTALL.MacOSX
index cbf8a014..8e36dd86 100644
--- a/INSTALL.MacOSX
+++ b/INSTALL.MacOSX
@@ -1,5 +1,33 @@
 COMPILING AND INSTALLING
 
+The following instructions apply for use of package installer `homebrew`. Legacy instructions below were also showing
+use of MacPorts for older version of python and it has not been tested
+
+# Installation with `homebrew`
+
+Install brew (https://docs.brew.sh/Installation) and Xcode ((from the App Store or from the apple developer website), agree to the Xcode license and install the Xcode command-line tools.
+
+Install the following brew package:
+$ brew install texlive
+$ brew install boost boost-python3
+$ brew install gcc lapack blitz scons
+
+Use pyenv (https://github.com/pyenv/pyenv) to install your favorite python 3 version.
+
+install CAMFR python requirements:
+
+$ pip install requirements.txt
+
+$ cp machine_cfg.py.MacOSX-brew machine_cfg.py
+
+check paths in machine_cfg.py - and run installation
+
+$ pip install .
+
+# Installation with MacPorts - and python 2.7
+
+(this won't work on Apple M1 where python 2.7 is deprecated)
+
 For the compilation of CAMFR on MacOSX, it is relatively straightforward to use MacPorts to install all dependencies, although with the drawback that the built-in Mac OS installation of Python will be superseded by the MacPorts python installation.
 The following instructions will install a new Python 2.7 interpreter, alongside your built-in MacOS Python installation. Be sure to used the correct python shell when running the CAMFR install script, and note that CAMFR will only be installed to this particular python interpreter unless you manually copy the module into another `site-packages` folder.
 
diff --git a/camfr/camfr_wrap.cpp b/camfr/camfr_wrap.cpp
index 37361684..10b66ee3 100644
--- a/camfr/camfr_wrap.cpp
+++ b/camfr/camfr_wrap.cpp
@@ -365,11 +365,12 @@ inline void free_tmp_interfaces(Waveguide& w)
 struct cVector_to_python
 {
   static PyObject* convert(const cVector& c)
-  {    
-    int dim[1]; dim[0] = c.rows();
+  {
+    npy_intp dim[1];
+    dim[0] = c.rows();
 
-    PyArrayObject* result 
-      = (PyArrayObject*) PyArray_FromDims(1, dim, PyArray_CDOUBLE);
+    PyArrayObject* result
+      = (PyArrayObject*) PyArray_SimpleNew(1, dim, PyArray_CDOUBLE);
 
     for (int i=0; i<c.rows(); i++)
         *(Complex*)(result->data + i*result->strides[0]) = c(i+1);
@@ -427,10 +428,10 @@ struct cMatrix_to_python
 {
   static PyObject* convert(const cMatrix& c)
   {
-    int dim[2]; dim[0] = c.rows(); dim[1] = c.columns();
+    npy_intp dim[2]; dim[0] = c.rows(); dim[1] = c.columns();
   
     PyArrayObject* result 
-      = (PyArrayObject*) PyArray_FromDims(2, dim, PyArray_CDOUBLE);
+      = (PyArrayObject*) PyArray_SimpleNew(2, dim, PyArray_CDOUBLE);
     
     for (int i=0; i<c.rows(); i++)
       for (int j=0; j<c.columns(); j++)
diff --git a/camfr/math/calculus/croot/allroots.cpp b/camfr/math/calculus/croot/allroots.cpp
index cb052dd4..f229c825 100755
--- a/camfr/math/calculus/croot/allroots.cpp
+++ b/camfr/math/calculus/croot/allroots.cpp
@@ -64,7 +64,7 @@ vector<Complex> roots_contour(const Contour& contour,
   // Check if all G's are zero.
 
   bool all_zeros = true;
-  for (unsigned int i=0; i<=2*N-1; i++)
+  for (unsigned int i=0; i<G.size(); i++)
     if (abs(G[i]) > 1e-10)
       all_zeros = false;
 
@@ -83,7 +83,10 @@ vector<Complex> roots_contour(const Contour& contour,
       A(i,j) = G[i+j-2];
 
   for (int i=1; i<=N; i++)
-    B(i,1) = -G[N+i-1];
+    // CHECK: something strange here, I added boundary checking - however logically we should start from where the
+    // previous loop ends with N+i-2 - but changing this does not work, so I am adding this boundary check
+    if (N+i-1!=G.size())
+        B(i,1) = -G[N+i-1];
 
   // Solve linear system and exploit the symmetry. Note: the fact that this
   // is a Hankel matrix is not used, since speedups would probably be minor
diff --git a/camfr/math/calculus/croot/patterson_z_n.cpp b/camfr/math/calculus/croot/patterson_z_n.cpp
index 04b2c2b1..0d06e71e 100644
--- a/camfr/math/calculus/croot/patterson_z_n.cpp
+++ b/camfr/math/calculus/croot/patterson_z_n.cpp
@@ -43,6 +43,7 @@ vector<Complex> patterson_z_n(ComplexFunction& f,
   {
     vector<Complex> result;
 
+    // TODO: - missing abs_error definition
     for (unsigned int i=0; i<=M; i++)
       result.push_back(0.0);
     
@@ -277,10 +278,10 @@ vector<Complex> patterson_quad_z_n_sub
   // Do adaptive subdivision of interval.
 
   vector<Complex> result1
-    = patterson_quad_z_n_sub(f, a, (a+b)/2., M,eps,mu,result_estimate,max_k);
+    = patterson_quad_z_n_sub(f, a, (a+b)/2., result.size()-1,eps,mu,result_estimate,max_k);
 
   vector<Complex> result2
-    = patterson_quad_z_n_sub(f, (a+b)/2., b, M,eps,mu,result_estimate,max_k);
+    = patterson_quad_z_n_sub(f, (a+b)/2., b, result.size()-1,eps,mu,result_estimate,max_k);
 
   unsigned int new_M
     = (result1.size()>result2.size()) ? result2.size() : result1.size();
@@ -314,12 +315,12 @@ vector<Complex> patterson_quad_z_n(ComplexFunction& f,
     return result;
 
   // Do adaptive subdivision of interval
-  
+
   vector<Complex> result1
-    = patterson_quad_z_n_sub(f, a, (a+b)/2., M, eps, mu, result, max_k);
+    = patterson_quad_z_n_sub(f, a, (a+b)/2., result.size()-1, eps, mu, result, max_k);
 
   vector<Complex> result2
-    = patterson_quad_z_n_sub(f, (a+b)/2., b, M, eps, mu, result, max_k);
+    = patterson_quad_z_n_sub(f, (a+b)/2., b, result.size()-1, eps, mu, result, max_k);
 
   unsigned int new_M
     = (result1.size()>result2.size()) ? result2.size() : result1.size();
diff --git a/camfr/math/calculus/polyroot/polyroot.cpp b/camfr/math/calculus/polyroot/polyroot.cpp
index fe01da01..a0e5ee52 100644
--- a/camfr/math/calculus/polyroot/polyroot.cpp
+++ b/camfr/math/calculus/polyroot/polyroot.cpp
@@ -55,7 +55,6 @@ vector<Complex> polyroot(const vector<Complex>& coef)
     py_error("Warning: polyroot solver did not converge.");
 
   delete [] coef_r; delete [] coef_i;
-  delete [] root_r, delete [] root_i;
 
   // Return results.
 
@@ -63,5 +62,7 @@ vector<Complex> polyroot(const vector<Complex>& coef)
   for (unsigned int i=0; i<N; i++)
     results.push_back(Complex(root_r[i],root_i[i]));
 
+  delete [] root_r, delete [] root_i;
+
   return results;
 }
diff --git a/camfr/mode.h b/camfr/mode.h
index 4dfb42c3..ed745cb5 100644
--- a/camfr/mode.h
+++ b/camfr/mode.h
@@ -171,11 +171,12 @@ struct modesorter_BDM
       const double ra2 = real(a->get_kz() * a->get_kz());
       const double rb2 = real(b->get_kz() * b->get_kz());
       
-      if ( (ra2 < 0) && (rb2 < 0) )
+      if ( (ra2 < 0) && (rb2 < 0) ) {
         if ( (a->pol != TE) && (a->pol != TM) )
           return ( ra2 < rb2 );
         else
           return ( ra2 > rb2 );
+      }
 
       if ( (ra2 > 0) && (rb2 > 0) )
         return ( ra2 > rb2 );
diff --git a/camfr/primitives/blochsection/blochsection.cpp b/camfr/primitives/blochsection/blochsection.cpp
index bc43c2a7..4145a33a 100644
--- a/camfr/primitives/blochsection/blochsection.cpp
+++ b/camfr/primitives/blochsection/blochsection.cpp
@@ -1287,6 +1287,8 @@ int UniformBlochSection::order(Polarisation pol, int Mx, int My) const
     if ( (m->pol == pol) && (m->get_Mx() == Mx) && (m->get_My() == My) )
       return i;
   }
+  // CHECK - what should be default return?
+  return 0;
 }
 
 
diff --git a/camfr/primitives/slab/isoslab/slab.cpp b/camfr/primitives/slab/isoslab/slab.cpp
index d982804b..38f854b5 100755
--- a/camfr/primitives/slab/isoslab/slab.cpp
+++ b/camfr/primitives/slab/isoslab/slab.cpp
@@ -1489,7 +1489,7 @@ void Slab_M::build_modeset(vector<Complex>& kt)
 
     if (i == 0)
     {
-      if (real(materials[i]->eps_mu()) > real(materials[i+1]->eps_mu()))
+      if (materials.size()>1 && real(materials[i]->eps_mu()) > real(materials[i+1]->eps_mu()))
         is_core = true;
     }
     else if (i == materials.size()-1)
diff --git a/camfr/stack.cpp b/camfr/stack.cpp
index 5c88e305..a83750e8 100644
--- a/camfr/stack.cpp
+++ b/camfr/stack.cpp
@@ -109,16 +109,18 @@ StackImpl::StackImpl(const Expression& e_, unsigned int no_of_periods_)
   }
   
   // Create chunks.
-  
+
   for (unsigned int i=0; i<e.get_size(); i++)
   {
     Term* t1 = e.get_term(i);
-    Term* t2 = e.get_term(i+1);
+    Term* t2 = 0;
+    if (i+1 < e.get_size())
+        t2 = e.get_term(i+1);
 
     // No waveguide.
     
-    if (t1->get_type() != WAVEGUIDE)
-      if ( (i+1 < e.get_size()) && (t2->get_type() == WAVEGUIDE) )
+    if (t1->get_type() != WAVEGUIDE) {
+      if (t2 && (t2->get_type() == WAVEGUIDE) )
       {
         Complex d = t2->get_d();
 
@@ -139,7 +141,8 @@ StackImpl::StackImpl(const Expression& e_, unsigned int no_of_periods_)
       }
       else
         chunks.push_back(Chunk(t1->get_sc(), 0.0));
-    
+    }
+
     // Waveguide.
 
     if (t1->get_type() == WAVEGUIDE)
diff --git a/camfr/util/cvector.cpp b/camfr/util/cvector.cpp
index b4d833a3..54dfe2cf 100755
--- a/camfr/util/cvector.cpp
+++ b/camfr/util/cvector.cpp
@@ -25,8 +25,13 @@ vector<Complex> operator+(const vector<Complex>& a,
 {  
   vector<Complex> result;
 
-  for (unsigned int i=0; i<a.size(); i++)
-    result.push_back(a[i]+b[i]);
+  for (unsigned int i=0; i<std::max(a.size(), b.size()); i++)
+    if (i<a.size() && i<b.size())
+        result.push_back(a[i]+b[i]);
+    else if (i<a.size())
+        result.push_back(a[i]);
+    else
+        result.push_back(b[i]);
 
   return result; 
 }
@@ -51,10 +56,14 @@ vector<Complex> operator-(const vector<Complex>& a,
 {  
   vector<Complex> result;
 
-  for (unsigned int i=0; i<a.size(); i++)
-    result.push_back(a[i]-b[i]);
-
-  return result; 
+  for (unsigned int i=0; i<std::max(a.size(), b.size()); i++)
+    if (i<a.size() && i<b.size())
+        result.push_back(a[i]-b[i]);
+    else if (i<a.size())
+        result.push_back(a[i]);
+    else
+        result.push_back(-b[i]);
+  return result;
 }
 
 vector<Complex>& operator-=(      vector<Complex>& a,
diff --git a/camfr/waveguide.h b/camfr/waveguide.h
index 45d2c859..aec52214 100644
--- a/camfr/waveguide.h
+++ b/camfr/waveguide.h
@@ -61,10 +61,12 @@ class Waveguide
 {
   public:
 
-    Waveguide(bool uniform_=false, Material* core_=NULL)
-      : uniform(uniform_), core(core_) {}
-    
-    virtual ~Waveguide() {interface_cache.deregister(this);} // Important!
+    Waveguide(bool uniform_=false, Material* core_=NULL): uniform(uniform_), core(new Material(*core_)) {}
+
+    virtual ~Waveguide() {
+        delete core;
+        interface_cache.deregister(this);  // Important!
+    }
 
     bool    is_uniform() const {return uniform;}
     Material* get_core() const {return core;}
diff --git a/machine_cfg.py.MacOSX-brew b/machine_cfg.py.MacOSX-brew
index 2dea44c2..9881211f 100644
--- a/machine_cfg.py.MacOSX-brew
+++ b/machine_cfg.py.MacOSX-brew
@@ -8,8 +8,14 @@
 # brew install boost-python3
 # brew install texlive
 
+# make sure to have tcl-tk installed (brew install tcl-tk) before installing python (important because CAMFR has some built-in dependencies with tkinter)
 # define python version with pyenv and define virtualenv with venv
 
+# finally update the variable PATHTOPYENVPYTHONVERSION/.pyenv..., with path to your pyenv python version (for Python.h include) and .venv version (for numpy includes) /PATHTOVENV/venv/lib/python3.9/...
+
+# check compilation with "make"
+# then you can install with pip install .
+
 #
 # is building on ARM mac book pro Monterey 12.1
 
@@ -45,7 +51,7 @@ flags = fflags + " -ftemplate-depth-60"
 include_dirs = ["/opt/local/include",
                 "/opt/homebrew/Cellar//boost/1.78.0_1/include/",
                 "/opt/homebrew/Cellar/blitz/1.0.2/include/",
-                "~/.pyenv/versions/3.9.11/include/python3.9", "venv/lib/python3.9/site-packages"]
+                "/PATHTOPYENVPYTHONVERSION/.pyenv/versions/3.9.11/include/python3.9", "/PATHTOVENV/venv/lib/python3.9/site-packages"]
 
 library_dirs = [ "/opt/homebrew/Cellar//boost/1.78.0_1/lib", "/opt/homebrew/Cellar/boost-python3/1.78.0/lib",
                  "/opt/homebrew/Cellar/blitz/1.0.2/lib","/opt//homebrew/Cellar/gcc/11.3.0/lib/gcc/11/"
diff --git a/requirements.txt b/requirements.txt
index ed2733f9..06a4c992 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -2,7 +2,6 @@ Image==1.5.33
 ImageFilter==0.1.0
 matplotlib==3.5.2
 MLab==1.1.4
-Numeric==24.2
 numpy==1.22.3
 Pillow==9.1.0
 pymat==0.0.2