diff --git a/inst/@sym/private/elementwise_op.m b/inst/@sym/private/elementwise_op.m
index fe68ca95e..d03e0e94a 100644
--- a/inst/@sym/private/elementwise_op.m
+++ b/inst/@sym/private/elementwise_op.m
@@ -1,4 +1,4 @@
-%% Copyright (C) 2014, 2016, 2018-2019, 2022 Colin B. Macdonald
+%% Copyright (C) 2014, 2016, 2018-2019, 2021-2022 Colin B. Macdonald
 %% Copyright (C) 2016 Lagu
 %%
 %% This file is part of OctSymPy.
@@ -84,7 +84,7 @@
           % Make sure all matrices in the input are the same size, and set q to one of them
           'q = None'
           'for A in _ins:'
-          '    if isinstance(A, MatrixBase):'
+          '    if isinstance(A, (MatrixBase, NDimArray)):'
           '        if q is None:'
           '            q = A'
           '        else:'
@@ -97,7 +97,10 @@
           'elements = []'
           'for i in range(0, len(q)):'
           '    elements.append(_op(*[k[i] if isinstance(k, MatrixBase) else k for k in _ins]))'
-          'return Matrix(*q.shape, elements)' ];
+          'if all(isinstance(x, Expr) for x in elements):'
+          '    return Matrix(*q.shape, elements)'
+          'dbout(f"elementwise_op: returning an Array not a Matrix")'
+          'return Array(elements, shape=q.shape)' ];
 
   z = pycall_sympy__ (cmd, varargin{:});
 
diff --git a/inst/private/python_header.py b/inst/private/python_header.py
index 925e7a50b..05d016c05 100644
--- a/inst/private/python_header.py
+++ b/inst/private/python_header.py
@@ -152,7 +152,7 @@ def octoutput(x, et):
             f.text = str(OCTCODE_BOOL)
             f = ET.SubElement(a, "f")
             f.text = str(x)
-        elif x is None or isinstance(x, (sp.Basic, sp.MatrixBase)):
+        elif x is None or isinstance(x, (sp.Basic, sp.MatrixBase, sp.NDimArray)):
             # FIXME: is it weird to pretend None is a SymPy object?
             if isinstance(x, (sp.Matrix, sp.ImmutableMatrix)):
                 _d = x.shape
@@ -161,6 +161,9 @@ def octoutput(x, et):
                 _d = [float(r) if (isinstance(r, sp.Basic) and r.is_Integer)
                       else float('nan') if isinstance(r, sp.Basic)
                       else r for r in x.shape]
+            elif isinstance(x, sp.NDimArray):
+                _d = x.shape
+                dbout(f"I am here with an array with shape {_d}")
             elif x is None:
                 _d = (1,1)
             else: