diff --git a/geopdes/inst/msh/@msh_cartesian/msh_eval_boundary_side.m b/geopdes/inst/msh/@msh_cartesian/msh_eval_boundary_side.m
index 67bb2e4c..03b6ee6e 100644
--- a/geopdes/inst/msh/@msh_cartesian/msh_eval_boundary_side.m
+++ b/geopdes/inst/msh/@msh_cartesian/msh_eval_boundary_side.m
@@ -154,7 +154,7 @@
     % normal direction (the boundary manifold is not enough).
     element_list = element_list(:)';
     indices = cell (msh_side.ndim, 1);
-    [indices{:}] = ind2sub (msh_side.nel_dir, element_list);
+    [indices{:}] = ind2sub ([msh_side.nel_dir, 1], element_list);
     indices = cell2mat (indices);
     
     qn_elems = arrayfun(@(ii) {msh.boundary(iside).qn{ii}(:,indices(ii,:))}, 1:msh_side.ndim);