diff --git a/demo/data/Head03_16kHz.ply b/demo/data/Head03_16kHz.ply
deleted file mode 100755
index db0e6e39416f5c1801e25e9924737235cc3ba8a9..0000000000000000000000000000000000000000
Binary files a/demo/data/Head03_16kHz.ply and /dev/null differ
diff --git a/demo/data/Head03_20kHz.ply b/demo/data/Head03_20kHz.ply
deleted file mode 100755
index 942a72f22a05d07a4cbe33b7b764751a4423ee6c..0000000000000000000000000000000000000000
Binary files a/demo/data/Head03_20kHz.ply and /dev/null differ
diff --git a/demo/dirichlet.cpp b/demo/dirichlet.cpp
index 5b119a8a5c781fbfbf4866b0ad277a8575cc496d..d48c252928b77fd9aae6df2b134be3ecd21c1e86 100644
--- a/demo/dirichlet.cpp
+++ b/demo/dirichlet.cpp
@@ -99,8 +99,8 @@ int main(int argc, char* argv[])
     {
         return singleLayer<std::complex<double>>(Pvtx,u,k,Ix,Iy);
     };
-    matrix<std::complex<double>> Pdom = Pinc;
-    hmatrix<std::complex<double>> Sdom(Pvtx,u.dof(),tol,Sdomfct,lambda,Pdom);
+    hmatrix<std::complex<double>> Sdom(Pvtx,u.dof(),tol,Sdomfct);
+    matrix<std::complex<double>> Pdom = mtimes(Sdom,lambda) + Pinc;
     toc();
     
     // Graphical representation
@@ -112,12 +112,12 @@ int main(int argc, char* argv[])
     tic();
     matrix<double> theta = transpose(linspace(0,2*M_PI,1e3));
     matrix<double> Nu = horzcat(horzcat(cos(theta),sin(theta)),zeros(size(theta)));
-    auto Sinffct = [&Nu,&u,&k](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
+    matrix<std::complex<double>> Pinf(numel(theta),1);
+    for (std::size_t i=0; i<numel(theta); ++i)
     {
-        return singleLayerInfinite<std::complex<double>>(Nu,u,k,Ix,Iy);
-    };
-    matrix<std::complex<double>> Pinf = zeros(size(Nu,1),size(lambda,2));
-    hmatrix<std::complex<double>> Sinf(Nu,u.dof(),0,Sinffct,lambda,Pinf);
+        auto Sinf = singleLayerInfinite<std::complex<double>>(Nu,u,k,i);
+        Pinf(i) = (std::complex<double>) mtimes(Sinf,lambda);
+    }
     toc();
 
     // Graphical representation
diff --git a/demo/dirichletBW.cpp b/demo/dirichletBW.cpp
index c95e45fb29d2c0cc59f21f3a3c63b0cfad24e85f..7ad77f0df79d2ef5e5795b444f56099485cd100f 100644
--- a/demo/dirichletBW.cpp
+++ b/demo/dirichletBW.cpp
@@ -82,13 +82,13 @@ int main (int argc, char* argv[])
     {
         return singleLayer<std::complex<double>>(v,u,k,Ix,Iy);
     };
+    hmatrix<std::complex<double>> Sbnd(v.dof(),u.dof(),tol,Sbndfct);
     auto Dbndfct = [&v,&u,&k,&Id](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
     {
         return 0.5*eval(Id(Ix,Iy)) + doubleLayer<std::complex<double>>(v,u,k,Ix,Iy);
     };
-    matrix<std::complex<double>> Pbnd = zeros(size(v.dof(),1),size(mu,2));
-    hmatrix<std::complex<double>> Sbnd(v.dof(),u.dof(),tol,Sbndfct,lambda,Pbnd);
-    hmatrix<std::complex<double>> Dbnd(v.dof(),u.dof(),tol,Dbndfct,-mu,Pbnd);
+    hmatrix<std::complex<double>> Dbnd(v.dof(),u.dof(),tol,Dbndfct);
+    matrix<std::complex<double>> Pbnd = mtimes(Sbnd,lambda) - mtimes(Dbnd,mu);
     toc();
     Pbnd = gmres(Id,Pbnd,tol,100) + planeWave(v.dof(),U,k);
     
@@ -98,13 +98,13 @@ int main (int argc, char* argv[])
     {
         return singleLayer<std::complex<double>>(Pvtx,u,k,Ix,Iy);
     };
+    hmatrix<std::complex<double>> Sdom(Pvtx,u.dof(),tol,Sdomfct);
     auto Ddomfct = [&Pvtx,&u,&k](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
     {
         return doubleLayer<std::complex<double>>(Pvtx,u,k,Ix,Iy);
     };
-    matrix<std::complex<double>> Pdom = Pinc;
-    hmatrix<std::complex<double>> Sdom(Pvtx,u.dof(),tol,Sdomfct,lambda,Pdom);
-    hmatrix<std::complex<double>> Ddom(Pvtx,u.dof(),tol,Ddomfct,-mu,Pdom);
+    hmatrix<std::complex<double>> Ddom(Pvtx,u.dof(),tol,Ddomfct);
+    matrix<std::complex<double>> Pdom = mtimes(Sdom,lambda) - mtimes(Ddom,mu) + Pinc;
     toc();
     for (std::size_t l=0; l<size(Pvtx,1); ++l)
     {
@@ -120,17 +120,13 @@ int main (int argc, char* argv[])
     tic();
     matrix<double> theta = transpose(linspace(0,2*M_PI,1e3));
     matrix<double> Nu = horzcat(horzcat(cos(theta),sin(theta)),zeros(size(theta)));
-    auto Sinffct = [&Nu,&u,&k](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
+    matrix<std::complex<double>> Pinf(numel(theta),1);
+    for (std::size_t i=0; i<numel(theta); ++i)
     {
-        return singleLayerInfinite<std::complex<double>>(Nu,u,k,Ix,Iy);
-    };
-    auto Dinffct = [&Nu,&u,&k](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
-    {
-        return doubleLayerInfinite<std::complex<double>>(Nu,u,k,Ix,Iy);
-    };
-    matrix<std::complex<double>> Pinf = zeros(size(Nu,1),size(lambda,2));
-    hmatrix<std::complex<double>> Sinf(Nu,u.dof(),0,Sinffct,lambda,Pinf);
-    hmatrix<std::complex<double>> Dinf(Nu,u.dof(),0,Dinffct,-mu,Pinf);
+        auto Sinf = singleLayerInfinite<std::complex<double>>(Nu,u,k,i);
+        auto Dinf = doubleLayerInfinite<std::complex<double>>(Nu,u,k,i);
+        Pinf(i) = (std::complex<double>) (mtimes(Sinf,lambda) - mtimes(Dinf,mu));
+    }
     toc();
     
     // Graphical representation
diff --git a/demo/head.cpp b/demo/head.cpp
index fe18bfa92c4a9d18669188241fe33d1ea8616bd7..15adf82e00cc9e4ad3669e9051404e4d1d648962 100644
--- a/demo/head.cpp
+++ b/demo/head.cpp
@@ -75,14 +75,14 @@ int main (int argc, char* argv[])
     {
         return 0.5*eval(Id(Ix,Iy)) + doubleLayer<std::complex<double>>(v,u,k,Ix,Iy);
     };
-    matrix<std::complex<double>> Pbnd = zeros(size(v.dof(),1),size(mu,2));
-    hmatrix<std::complex<double>> Dbnd(v.dof(),u.dof(),tol,Dbndfct,mu,Pbnd);
+    hmatrix<std::complex<double>> Dbnd(v.dof(),u.dof(),tol,Dbndfct);
+    matrix<std::complex<double>> Pbnd = mtimes(Dbnd,mu);
     toc();
     Pbnd = - gmres(Id,Pbnd,tol,100) + planeWave(v.dof(),U,k);
     
     // Graphical representation
     figure fig2;
-    trimesh(fig2,Stri,Svtx,real(Pbnd));
+    trimesh(fig2,Stri,Svtx,abs(Pbnd));
     
     // Export in .vtk file
     triwrite("/Users/aussal/Desktop/resultats/","head.vtk",Stri,Svtx,real(Pbnd));
diff --git a/demo/neumann.cpp b/demo/neumann.cpp
index b11ab936e9587b717c583a5e5e3ef4e2c90ed68e..6c41a00f6fdc4472f1ee6e0c2341dfd954d27fef 100644
--- a/demo/neumann.cpp
+++ b/demo/neumann.cpp
@@ -83,8 +83,8 @@ int main (int argc, char* argv[])
     {
         return 0.5*eval(Id(Ix,Iy)) + doubleLayer<std::complex<double>>(v,u,k,Ix,Iy);
     };
-    matrix<std::complex<double>> Pbnd = zeros(size(v.dof(),1),size(mu,2));
-    hmatrix<std::complex<double>> Dbnd(v.dof(),u.dof(),tol,Dbndfct,mu,Pbnd);
+    hmatrix<std::complex<double>> Dbnd(v.dof(),u.dof(),tol,Dbndfct);
+    matrix<std::complex<double>> Pbnd = mtimes(Dbnd,mu);
     toc();
     Pbnd = - gmres(Id,Pbnd,tol,100) + planeWave(v.dof(),U,k);
 
@@ -94,8 +94,9 @@ int main (int argc, char* argv[])
     {
         return doubleLayer<std::complex<double>>(Pvtx,u,k,Ix,Iy);
     };
-    matrix<std::complex<double>> Pdom = Pinc;
-    hmatrix<std::complex<double>> Ddom(Pvtx,u.dof(),tol,Ddomfct,-mu,Pdom);
+    hmatrix<std::complex<double>> Ddom(Pvtx,u.dof(),tol,Ddomfct);
+    matrix<std::complex<double>> Pdom = - mtimes(Ddom,mu) + Pinc;
+    
     toc();
     for (std::size_t l=0; l<size(Pvtx,1); ++l)
     {
@@ -111,12 +112,12 @@ int main (int argc, char* argv[])
     tic();
     matrix<double> theta = transpose(linspace(0,2*M_PI,1e3));
     matrix<double> Nu = horzcat(horzcat(cos(theta),sin(theta)),zeros(size(theta)));
-    auto Dinffct = [&Nu,&u,&k](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
+    matrix<std::complex<double>> Pinf(numel(theta),1);
+    for (std::size_t i=0; i<numel(theta); ++i)
     {
-        return doubleLayerInfinite<std::complex<double>>(Nu,u,k,Ix,Iy);
-    };
-    matrix<std::complex<double>> Pinf = zeros(size(Nu,1),size(mu,2));
-    hmatrix<std::complex<double>> Dinf(Nu,u.dof(),0,Dinffct,-mu,Pinf);
+        auto Dinf = doubleLayerInfinite<std::complex<double>>(Nu,u,k,i);
+        Pinf(i) = (std::complex<double>) -mtimes(Dinf,mu);
+    }
     toc();
     
     // Graphical representation
diff --git a/demo/neumannBW.cpp b/demo/neumannBW.cpp
index ac94c7774cbb95e28216e22cce4e7eb3b928fad5..7d992ea11f917891fc9702f70ab69d6e7279142e 100644
--- a/demo/neumannBW.cpp
+++ b/demo/neumannBW.cpp
@@ -82,13 +82,13 @@ int main (int argc, char* argv[])
     {
         return singleLayer<std::complex<double>>(v,u,k,Ix,Iy);
     };
+    hmatrix<std::complex<double>> Sbnd(v.dof(),u.dof(),tol,Sbndfct);
     auto Dbndfct = [&v,&u,&k,&Id](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
     {
         return 0.5*eval(Id(Ix,Iy)) + doubleLayer<std::complex<double>>(v,u,k,Ix,Iy);
     };
-    matrix<std::complex<double>> Pbnd = zeros(size(v.dof(),1),size(mu,2));
-    hmatrix<std::complex<double>> Sbnd(v.dof(),u.dof(),tol,Sbndfct,lambda,Pbnd);
-    hmatrix<std::complex<double>> Dbnd(v.dof(),u.dof(),tol,Dbndfct,-mu,Pbnd);
+    hmatrix<std::complex<double>> Dbnd(v.dof(),u.dof(),tol,Dbndfct);
+    matrix<std::complex<double>> Pbnd = mtimes(Sbnd,lambda) - mtimes(Dbnd,mu);
     toc();
     Pbnd = gmres(Id,Pbnd,tol,100) + planeWave(v.dof(),U,k);
     
@@ -98,13 +98,13 @@ int main (int argc, char* argv[])
     {
         return singleLayer<std::complex<double>>(Pvtx,u,k,Ix,Iy);
     };
+    hmatrix<std::complex<double>> Sdom(Pvtx,u.dof(),tol,Sdomfct);
     auto Ddomfct = [&Pvtx,&u,&k](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
     {
         return doubleLayer<std::complex<double>>(Pvtx,u,k,Ix,Iy);
     };
-    matrix<std::complex<double>> Pdom = Pinc;
-    hmatrix<std::complex<double>> Sdom(Pvtx,u.dof(),tol,Sdomfct,lambda,Pdom);
-    hmatrix<std::complex<double>> Ddom(Pvtx,u.dof(),tol,Ddomfct,-mu,Pdom);
+    hmatrix<std::complex<double>> Ddom(Pvtx,u.dof(),tol,Ddomfct);
+    matrix<std::complex<double>> Pdom = mtimes(Sdom,lambda) - mtimes(Ddom,mu) + Pinc;
     toc();
     for (std::size_t l=0; l<size(Pvtx,1); ++l)
     {
@@ -120,17 +120,13 @@ int main (int argc, char* argv[])
     tic();
     matrix<double> theta = transpose(linspace(0,2*M_PI,1e3));
     matrix<double> Nu = horzcat(horzcat(cos(theta),sin(theta)),zeros(size(theta)));
-    auto Sinffct = [&Nu,&u,&k](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
+    matrix<std::complex<double>> Pinf(numel(theta),1);
+    for (std::size_t i=0; i<numel(theta); ++i)
     {
-        return singleLayerInfinite<std::complex<double>>(Nu,u,k,Ix,Iy);
-    };
-    auto Dinffct = [&Nu,&u,&k](matrix<std::size_t> Ix, matrix<std::size_t> Iy)
-    {
-        return doubleLayerInfinite<std::complex<double>>(Nu,u,k,Ix,Iy);
-    };
-    matrix<std::complex<double>> Pinf = zeros(size(Nu,1),size(lambda,2));
-    hmatrix<std::complex<double>> Sinf(Nu,u.dof(),0,Sinffct,lambda,Pinf);
-    hmatrix<std::complex<double>> Dinf(Nu,u.dof(),0,Dinffct,-mu,Pinf);
+        auto Sinf = singleLayerInfinite<std::complex<double>>(Nu,u,k,i);
+        auto Dinf = doubleLayerInfinite<std::complex<double>>(Nu,u,k,i);
+        Pinf(i) = (std::complex<double>) (mtimes(Sinf,lambda) - mtimes(Dinf,mu));
+    }
     toc();
     
     // Graphical representation