diff --git a/sympde/calculus/__init__.py b/sympde/calculus/__init__.py
index d468d7df..1a6a6448 100644
--- a/sympde/calculus/__init__.py
+++ b/sympde/calculus/__init__.py
@@ -1,3 +1,2 @@
 from .core     import *
 from .errors   import *
-from .matrices import *
diff --git a/sympde/calculus/core.py b/sympde/calculus/core.py
index ac59ef41..151c0c01 100644
--- a/sympde/calculus/core.py
+++ b/sympde/calculus/core.py
@@ -428,6 +428,14 @@ def __new__(cls, arg1, arg2, **options):
             a,b = b,a
             c   = -c
 
+        # after simplifications, the expression can be reduced to 0 in some cases
+        # here's an example;
+        # if b is a Vector and u and v are scalar functions,
+        # we should have the following
+        # cross(b*u,b*v) = u*v*cross(b,b) = 0
+        if a == b:
+            return S.Zero
+
         obj = Basic.__new__(cls, a, b)
 
         return c*obj
diff --git a/sympde/core/__init__.py b/sympde/core/__init__.py
index 4bc8bacf..8d4e49aa 100644
--- a/sympde/core/__init__.py
+++ b/sympde/core/__init__.py
@@ -1,3 +1,4 @@
 from .algebra      import *
 from .basic        import *
 from .utils        import *
+from .matrices     import *
diff --git a/sympde/calculus/matrices.py b/sympde/core/matrices.py
similarity index 74%
rename from sympde/calculus/matrices.py
rename to sympde/core/matrices.py
index 16b827e1..4cfe7c30 100644
--- a/sympde/calculus/matrices.py
+++ b/sympde/core/matrices.py
@@ -1,9 +1,15 @@
 from sympy                 import Expr, S
 from sympy                 import Add, Mul, Pow
+from sympy                 import Tuple
 from sympy                 import sympify
 from sympy.core.decorators import call_highest_priority
 from sympde.core.basic     import _coeffs_registery, Basic
 
+from sympy import ImmutableDenseMatrix
+
+import functools
+
+
 class MatrixSymbolicExpr(Expr):
     is_commutative = False
     _op_priority   = 20.0
@@ -289,6 +295,106 @@ def _sympystr(self, printer):
         sstr = printer.doprint
         return '{}[{}]'.format(sstr(self.args[0]),sstr(self.args[1]))
 
+class Matrix(MatrixSymbolicExpr):
+    def __new__(cls, mat, *, name):
+        if not isinstance(mat, list):
+            raise TypeError('Positional argument `mat` should be a list of lists.')
+
+        for row in mat:
+            if not isinstance(row, list):
+                raise TypeError('Each row of `mat` should be a list.')
+
+        row_sizes = {len(row) for row in mat}
+        if len(row_sizes) != 1:
+            raise ValueError('Each row of `mat` should have the same length.')
+
+        if not isinstance(name, str):
+            raise TypeError('Keyword-only argument `name` must be a string.')
+
+        # Call constructor of base class Expr with immutable arguments
+        new_mat = Tuple(*[Tuple(*row) for row in mat])
+        obj = Expr.__new__(cls, new_mat, name)
+        return obj
+
+    @property
+    def name(self):
+        return self.args[1]
+
+    def __getitem__(self, key):
+        i,j = key
+        return self.args[0][i][j]
+
+    def _sympystr(self, printer):
+        sstr = printer.doprint
+        return '{}'.format(sstr(self.name))
+
+    def __hash__(self):
+        return hash((self.name, self.args))
+
+    def __eq__(self, other):
+        if not isinstance(other, Matrix):
+            return False
+
+        if not self.name == other.name:
+            return False
+
+        result = all(x == y for x, y in zip(self.args[0], other.args[0]))
+        return result
+
+    def to_sympy(self):
+        return ImmutableDenseMatrix(self.args[0])
+
+
+class Vector(MatrixSymbolicExpr):
+
+    def __new__(cls, vec, *, name):
+        if not isinstance(vec, list):
+            raise TypeError('Positional argument `vec` should be a list.')
+
+        if not isinstance(name, str):
+            raise TypeError('Keyword-only argument `name` must be a string.')
+
+        # Call constructor of base class Expr with immutable arguments
+        new_vec = Tuple(*vec)
+        obj = Expr.__new__(cls, new_vec, name)
+        return obj
+
+    @property
+    def name(self):
+        return self.args[1]
+
+    def __getitem__(self, key):
+        i = key
+        return self.args[0][i]
+
+    def _sympystr(self, printer):
+        sstr = printer.doprint
+        return '{}'.format(sstr(self.name))
+
+    def __hash__(self):
+        return hash((self.name, self.args))
+
+    def __eq__(self, other):
+        # TODO BUG
+        # We should check that other is a Vector
+        # right now, there is a problem when using Cross
+        # see the linearity test of
+        # f = lambda u,v: cross(b*u, b*v)
+        # where b is a Vector
+        if not isinstance(other, Vector):
+            return False
+
+        if not self.name == other.name:
+            return False
+
+        result = self.args[0] == other.args[0]
+        return result
+
+    def to_sympy(self):
+        args = [[a] for a in self.args[0]]
+        return ImmutableDenseMatrix(args)
+
+
 Basic._constructor_postprocessor_mapping[MatrixSymbolicExpr] = {
     "Mul": [lambda x: MatSymbolicMul(*x.args)],
     "Add": [lambda x: MatSymbolicAdd(*x.args)]
diff --git a/sympde/expr/evaluation.py b/sympde/expr/evaluation.py
index e75189be..381c027e 100644
--- a/sympde/expr/evaluation.py
+++ b/sympde/expr/evaluation.py
@@ -20,12 +20,14 @@
                                  Dot_2d, Inner_2d, Cross_2d,
                                  Dot_3d, Inner_3d, Cross_3d)
 from sympde.core.utils import random_string
+from sympde.core.matrices import SymbolicDeterminant, Inverse, Transpose
+from sympde.core.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace
+from sympde.core.matrices import Matrix as SympdeMatrix
+from sympde.core.matrices import Vector as SympdeVector
 
 from sympde.calculus import jump, avg, minus, plus
 from sympde.calculus import Jump, is_zero
 from sympde.calculus.core import _generic_ops, _diff_ops
-from sympde.calculus.matrices import SymbolicDeterminant, Inverse, Transpose
-from sympde.calculus.matrices import MatSymbolicPow, MatrixElement, SymbolicTrace
 
 from sympde.topology.basic   import BasicDomain, Union, Interval
 from sympde.topology.basic   import Boundary, Interface
@@ -614,6 +616,12 @@ def eval(cls, expr, domain):
         elif isinstance(expr, Inverse):
             return cls.eval(expr.arg, domain=domain).inv()
 
+        elif isinstance(expr, SympdeMatrix):
+            return expr.to_sympy()
+
+        elif isinstance(expr, SympdeVector):
+            return expr.to_sympy()
+
         elif isinstance(expr, ScalarFunction):
             return expr
 
diff --git a/sympde/expr/expr.py b/sympde/expr/expr.py
index bd87dfba..c866fc8a 100644
--- a/sympde/expr/expr.py
+++ b/sympde/expr/expr.py
@@ -22,6 +22,8 @@
 from sympde.topology.space import ScalarFunction
 from sympde.topology.space import VectorFunction
 from sympde.topology.space import Trace, trace_0, trace_1
+from sympde.core.matrices import Matrix as SympdeMatrix
+from sympde.core.matrices import Vector
 
 from .errors import UnconsistentLinearExpressionError
 from .basic  import BasicForm
@@ -538,9 +540,12 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options):
         # ...
 
         # ...
-        is_vector = isinstance(expr, (Matrix, ImmutableDenseMatrix, Tuple, list, tuple))
+        is_vector = isinstance(expr, (Vector, SympdeMatrix, Matrix, ImmutableDenseMatrix, Tuple, list, tuple))
         if is_vector:
-            expr = ImmutableDenseMatrix(expr)
+            if isinstance(expr, (Vector, SympdeMatrix)):
+                expr = expr.to_sympy()
+            else:
+                expr = ImmutableDenseMatrix(expr)
         # ...
 
         # ...
@@ -610,9 +615,12 @@ def __new__(cls, expr, domain, kind='l2', evaluate=True, **options):
         # ...
 
         # ...
-        is_vector = isinstance(expr, (Matrix, ImmutableDenseMatrix, Tuple, list, tuple))
+        is_vector = isinstance(expr, (Vector, SympdeMatrix, Matrix, ImmutableDenseMatrix, Tuple, list, tuple))
         if is_vector:
-            expr = ImmutableDenseMatrix(expr)
+            if isinstance(expr, (Vector, SympdeMatrix)):
+                expr = expr.to_sympy()
+            else:
+                expr = ImmutableDenseMatrix(expr)
         # ...
 
         # ...
diff --git a/sympde/expr/tests/test_equation_2d.py b/sympde/expr/tests/test_equation_2d.py
index cb446f2f..5d66e10e 100644
--- a/sympde/expr/tests/test_equation_2d.py
+++ b/sympde/expr/tests/test_equation_2d.py
@@ -2,7 +2,6 @@
 
 #import pytest
 
-from sympy.core.containers import Tuple
 from sympy import pi, cos, sin
 from sympy import ImmutableDenseMatrix as Matrix
 
@@ -48,7 +47,7 @@ def test_equation_2d_1():
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
     int_2 = lambda expr: integral(B2, expr)
-    
+
     # ... bilinear/linear forms
     expr = dot(grad(v), grad(u))
     a1 = BilinearForm((v,u), int_0(expr))
@@ -169,7 +168,7 @@ def test_equation_2d_2():
     alpha = Constant('alpha', real=True)
 
     int_0 = lambda expr: integral(domain , expr)
-    
+
     s  = BilinearForm((tau,sigma), int_0(dot(grad(tau), grad(sigma))))
     m  = BilinearForm((tau,sigma), int_0(tau*sigma))
     b1 = BilinearForm((tau,dw), int_0(bracket(pn, dw) * tau))
@@ -202,7 +201,7 @@ def test_equation_2d_3():
     x,y = domain.coordinates
 
     B1 = Boundary(r'\Gamma_1', domain)
-    
+
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
 
@@ -235,7 +234,7 @@ def test_equation_2d_4():
     x,y = domain.coordinates
 
     B1 = Boundary(r'\Gamma_1', domain)
-    
+
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
 
@@ -283,7 +282,7 @@ def test_equation_2d_5():
     p,q = [element_of(V, name=i) for i in ['p', 'q']]
 
     int_0 = lambda expr: integral(domain , expr)
-    
+
     a0 = BilinearForm((v,u), int_0(inner(grad(v), grad(u))))
     print('     a0 done.')
     a1 = BilinearForm((q,p), int_0(p*q))
@@ -335,7 +334,7 @@ def test_equation_2d_6():
     u,v = [element_of(V, name=i) for i in ['u', 'v']]
 
     int_0 = lambda expr: integral(domain , expr)
-    
+
     # ...
     expr = kappa * dot(grad(u), grad(v)) + dot(b, grad(u)) * v
     a = BilinearForm((v,u), int_0(expr))
diff --git a/sympde/expr/tests/test_expr.py b/sympde/expr/tests/test_expr.py
index efca3a98..0231b6a4 100644
--- a/sympde/expr/tests/test_expr.py
+++ b/sympde/expr/tests/test_expr.py
@@ -4,15 +4,19 @@
 
 import pytest
 
-from sympy.core.containers import Tuple
 from sympy import Function
 from sympy import pi, cos, sin, exp
 from sympy import ImmutableDenseMatrix as Matrix
 
 from sympde.core     import Constant
-from sympde.calculus import grad, dot, inner, rot, div
+from sympde.calculus import dot, inner, outer, cross
+from sympde.calculus import grad, rot, div, curl
 from sympde.calculus import laplace, bracket, convect
 from sympde.calculus import jump, avg, Dn, minus, plus
+from sympde.core import Matrix as SympdeMatrix
+from sympde.core import Vector
+from sympde.core import Transpose
+from sympde.core import Vector
 
 from sympde.topology import dx1, dx2, dx3
 from sympde.topology import dx, dy, dz
@@ -33,6 +37,8 @@
 from sympde.expr.expr import Functional, Norm
 from sympde.expr.expr import linearize
 from sympde.expr.evaluation import TerminalExpr
+from sympde.expr.expr import is_linear_expression
+
 
 #==============================================================================
 def test_linear_expr_2d_1():
@@ -82,7 +88,7 @@ def test_linear_expr_2d_2():
     u,u1,u2 = [element_of(V, name=i) for i in ['u', 'u1', 'u2']]
     v,v1,v2 = [element_of(V, name=i) for i in ['v', 'v1', 'v2']]
 
-    g = Tuple(x,y)
+    g = Vector([x,y], name='g')
     l = LinearExpr(v, dot(g, v))
     print(l)
     print(l.expr)
@@ -93,8 +99,8 @@ def test_linear_expr_2d_2():
     # ...
 
     # ...
-    g1 = Tuple(x,0)
-    g2 = Tuple(0,y)
+    g1 = Vector([x,0], name='g1')
+    g2 = Vector([0,y], name='g2')
     l = LinearExpr((v1,v2), dot(g1, v1) + dot(g2, v2))
     print(l)
     print(l.expr)
@@ -132,7 +138,7 @@ def test_linear_form_2d_1():
     # ...
 
     # ...
-    g  = Tuple(x**2, y**2)
+    g = Vector([x**2, y**2], name='g')
     l  = LinearForm(v, int_1(v*dot(g, nn)))
     print(l)
 
@@ -142,7 +148,7 @@ def test_linear_form_2d_1():
     # ...
 
     # ...
-    g = Tuple(x**2, y**2)
+    g = Vector([x**2, y**2], name='g')
     l = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v))
 
     assert(len(l.domain.args) == 2)
@@ -162,7 +168,7 @@ def test_linear_form_2d_1():
     # ...
 
     # ...
-    g = Tuple(x,y)
+    g = Vector([x,y], name='g')
     l1 = LinearForm(u1, int_0(x*y*u1))
     l2 = LinearForm(u2, int_0(dot(grad(u2), g)))
 
@@ -213,7 +219,7 @@ def test_linear_form_2d_2():
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
 
-    g = Matrix((x,y))
+    g = Vector([x,y], name='g')
     l = LinearForm(v, int_0(dot(g, v)))
 
     assert(l.domain == domain.interior)
@@ -221,7 +227,7 @@ def test_linear_form_2d_2():
     # ...
 
     # ...
-    g = Matrix((x,y))
+    g = Vector([x,y], name='g')
     l1 = LinearForm(v1, int_0(dot(g, v1)))
     l = LinearForm(v, l1(v))
 
@@ -230,8 +236,8 @@ def test_linear_form_2d_2():
     # ...
 
     # ...
-    g1 = Matrix((x,0))
-    g2 = Matrix((0,y))
+    g1 = Vector([x,0], name='g1')
+    g2 = Vector([0,y], name='g2')
     l1 = LinearForm(v1, int_0(dot(v1, g1)))
     l2 = LinearForm(v2, int_0(dot(v2, g2)))
 
@@ -396,14 +402,14 @@ def test_terminal_expr_linear_2d_1():
     # ...
 
     # ...
-    g = Matrix((x**2, y**2))
+    g = Vector([x**2, y**2], name='g')
     l = LinearForm(v, int_1(v*dot(g, nn)))
     print(TerminalExpr(l, domain))
     print('')
     # ...
 
     # ...
-    g = Matrix((x**2, y**2))
+    g = Vector([x**2, y**2], name='g')
     l = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v))
     print(TerminalExpr(l, domain))
     print('')
@@ -417,7 +423,7 @@ def test_terminal_expr_linear_2d_1():
     # ...
 
     # ...
-    g = Matrix((x,y))
+    g = Vector([x,y], name='g')
     l1 = LinearForm(v1, int_0(x*y*v1))
     l2 = LinearForm(v2, int_0(dot(grad(v2), g)))
 
@@ -435,7 +441,7 @@ def test_terminal_expr_linear_2d_1():
     # ...
 
     # ...
-    g = Matrix((x**2, y**2))
+    g = Vector([x**2, y**2], name='g')
     l1 = LinearForm(v1, int_0(x*y*v1))
     l2 = LinearForm(v1, int_0(v1))
     l3 = LinearForm(v, int_1(v*dot(g, nn)))
@@ -464,14 +470,14 @@ def test_terminal_expr_linear_2d_2():
     int_0 = lambda expr: integral(domain , expr)
     int_1 = lambda expr: integral(B1, expr)
 
-    g = Matrix((x,y))
+    g = Vector([x, y], name='g')
     l = LinearForm(v, int_0(dot(g, v)))
     print(TerminalExpr(l, domain))
     print('')
     # ...
 
     # ...
-    g = Matrix((x,y))
+    g = Vector([x, y], name='g')
     l = LinearForm(v, int_0(dot(g, v) + div(v)))
     print(TerminalExpr(l, domain))
     print('')
@@ -479,14 +485,14 @@ def test_terminal_expr_linear_2d_2():
 
     # TODO
 #    # ...
-#    g = Tuple(x**2, y**2)
+#    g = Vector([x**2,y**2], name='g')
 #    l = LinearForm(v, v*trace_1(g, B1))
 #    print(TerminalExpr(l))
 #    print('')
 #    # ...
 #
 #    # ...
-#    g = Tuple(x**2, y**2)
+#    g = Vector([x**2,y**2], name='g')
 #    l = LinearForm(v, v*trace_1(g, B1) + x*y*v)
 #    print(TerminalExpr(l))
 #    print('')
@@ -500,7 +506,7 @@ def test_terminal_expr_linear_2d_2():
 #    # ...
 #
 #    # ...
-#    g = Tuple(x,y)
+#    g = Vector([x,y], name='g')
 #    l1 = LinearForm(v1, x*y*v1)
 #    l2 = LinearForm(v2, dot(grad(v2), g))
 #
@@ -518,7 +524,7 @@ def test_terminal_expr_linear_2d_2():
 #    # ...
 #
 #    # ...
-#    g = Tuple(x**2, y**2)
+#    g = Vector([x**2,y**2], name='g')
 #    l1 = LinearForm(v1, x*y*v1)
 #    l2 = LinearForm(v1, v1)
 #    l3 = LinearForm(v, v*trace_1(g, B1))
@@ -585,7 +591,7 @@ def test_terminal_expr_linear_2d_5():
 
     V      = ScalarFunctionSpace('V', domain)
 
-    B_neumann = [domain.get_boundary(axis=0, ext=-1), 
+    B_neumann = [domain.get_boundary(axis=0, ext=-1),
                 domain.get_boundary(axis=1, ext=-1)]
 
     if len(B_neumann) == 1:
@@ -1178,9 +1184,7 @@ def test_stabilization_2d_1():
     kappa = Constant('kappa', is_real=True)
     mu    = Constant('mu'   , is_real=True)
 
-    b1 = 1.
-    b2 = 0.
-    b = Matrix((b1, b2))
+    b = Vector([1., 0.], name='b')
 
     # right hand side
     f = x*y
@@ -1403,8 +1407,8 @@ def test_norm_2d_2():
     F = element_of(V, 'F')
 
     # ...
-    f = Matrix((sin(pi*x)*sin(pi*y), sin(pi*x)*sin(pi*y)))
-    expr = Matrix([F[0]-f[0], F[1]-f[1]])
+    f = Vector([sin(pi*x)*sin(pi*y), sin(pi*x)*sin(pi*y)], name='f')
+    expr = Vector([F[0]-f[0], F[1]-f[1]], name='expr')
     l2_norm_u = Norm(expr, domain, kind='l2')
     h1_norm_u = Norm(expr, domain, kind='h1')
 
@@ -1502,16 +1506,16 @@ def test_linearity_linear_form_2d_1():
     _ = LinearForm(v, int_0(x * y * v))
     _ = LinearForm(v, int_0(x * y * v + v))
 
-    g = Matrix((x**2, y**2))
+    g = Vector([x**2, y**2], name='g')
     _ = LinearForm(v, int_0(v * dot(g, nn)))
 
-    g = Matrix((x**2, y**2))
+    g = Vector([x**2, y**2], name='g')
     _ = LinearForm(v, int_1(v*dot(g, nn)) + int_0(x*y*v))
 
     l1 = LinearForm(v1, int_0(x * y * v1))
     _  = LinearForm(v, l1(v))
 
-    g  = Matrix((x,y))
+    g = Vector([x, y], name='g')
     l1 = LinearForm(v1, int_0(x*y*v1))
     l2 = LinearForm(v2, int_0(dot(grad(v2), g)))
     _  = LinearForm(v, l1(v) + l2(v))
@@ -1520,7 +1524,7 @@ def test_linearity_linear_form_2d_1():
     l2 = LinearForm(v1, int_0(v1))
     _  = LinearForm(v, l1(v) + kappa*l2(v))
 
-    g  = Matrix((x**2, y**2))
+    g = Vector([x**2, y**2], name='g')
     l1 = LinearForm(v1, int_0(x*y*v1))
     l2 = LinearForm(v1, int_0(v1))
     l3 = LinearForm(v, int_1(v*dot(g, nn)))
@@ -1878,8 +1882,8 @@ def test_interface_integral_4():
 
     assert len(expr) == 1
     assert expr[0].target == B.get_boundary(axis=0, ext=-1)
-    assert expr[0].expr[0,0]   == u2[0]*v2[0]
-    assert expr[0].expr[1,1]   == u2[1]*v2[1]
+    assert expr[0].expr[0,0] == u2[0]*v2[0]
+    assert expr[0].expr[1,1] == u2[1]*v2[1]
 
     a2 = BilinearForm((u2,v2), integral(I, dot(minus(u2),minus(v2))))
 
@@ -1887,10 +1891,223 @@ def test_interface_integral_4():
 
     assert len(expr) == 1
     assert expr[0].target == A.get_boundary(axis=0, ext=1)
-    assert expr[0].expr[0,0]   == u2[0]*v2[0]
-    assert expr[0].expr[1,1]   == u2[1]*v2[1]
+    assert expr[0].expr[0,0] == u2[0]*v2[0]
+    assert expr[0].expr[1,1] == u2[1]*v2[1]
+    # ...
+
+#==============================================================================
+def test_matrices_vectors():
+
+    domain = Domain('Omega', dim=3)
+    x,y,z = domain.coordinates
+
+    kappa = Constant('kappa', is_real=True)
+    mu    = Constant('mu',    is_real=True)
+    theta = Constant('theta', is_real=True)
+    a1    = Constant('a1',    is_real=True)
+    a2    = Constant('a2',    is_real=True)
+    b1    = Constant('b1',    is_real=True)
+    b2    = Constant('b2',    is_real=True)
+
+    V = ScalarFunctionSpace('V', domain)
+    W = VectorFunctionSpace('W', domain)
+
+    u, v  = elements_of(V, names='u, v')
+    F, G  = elements_of(W, names='F, G')
+
+    #A = SympdeMatrix([[cos(theta), 0], [0, sin(theta)]], name='A')
+    A = SympdeMatrix([[1, 2, 3], [3, 4, 5], [5, 6, 7]], name='A')
+    b = Vector([1, 2, 3], name='b')
+
+    # ...
+    b1 = b
+    assert b1 == b
+
+    b1 = Vector([1, 2, 3], name='b1')
+    assert not b1 == b
+    # ...
+
+    # ...
+    A1 = A
+    assert A1 == A
+
+    A1 = SympdeMatrix([[1, 2, 3], [3, 4, 5], [5, 6, 7]], name='A1')
+    assert not A1 == A
+    # ...
+
+    # ... some tests on cross operator & vectors
+    assert cross(b,b)     == 0
+    assert cross(b*u,b*u) == 0
+    assert cross(b*u,b*v) == 0
+    # ...
+
+    ########################################################################
+    #                    SCALAR-FUNCTION CASE
+    ########################################################################
+    # ...
+    f = lambda a: A*grad(u)
+
+    assert is_linear_expression(f(u), (u,))
+    # ...
 
     # ...
+    f = lambda u: b*u
+
+    assert is_linear_expression(f(u), (u,))
+    # ...
+
+    # ...
+    f = lambda u: dot(b,grad(u))
+
+    assert is_linear_expression(f(u), (u,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(A*grad(u), grad(v))
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(grad(u), A*grad(v))
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(A*grad(u), A*grad(v))
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(A.T*grad(u), grad(v))
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(grad(u), A.T*grad(v))
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(A.T*grad(u), A.T*grad(v))
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(A*grad(u), grad(v))
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(b,grad(u))*v
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(A*grad(u), grad(v)) + dot(b,grad(u))*v + dot(b,grad(v))*u
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: cross(b*u, b*v)
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: cross(b*u, grad(v))
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: cross(grad(u), b*v)
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda u,v: dot(cross(grad(u), b), cross(grad(v), b))
+
+    assert is_linear_expression(f(u,v), (u,))
+    assert is_linear_expression(f(u,v), (v,))
+    # ...
+
+    # ...
+    f = lambda F,G: dot(curl(cross(b, F)), cross(b, G))
+
+    assert is_linear_expression(f(F,G), (F,))
+    assert is_linear_expression(f(F,G), (G,))
+    # ...
+
+    ## ...
+    #f = lambda u,v: 0
+    #
+    #assert is_linear_expression(f(u,v), (u,))
+    #assert is_linear_expression(f(u,v), (v,))
+    ## ...
+    #
+    ## ...
+    #f = lambda u,v: 0
+    #
+    #assert is_linear_expression(f(u,v), (u,))
+    #assert is_linear_expression(f(u,v), (v,))
+    ## ...
+
+    ########################################################################
+    #                    VECTOR-FUNCTION CASE
+    ########################################################################
+    # ...
+    f = lambda F,G: cross(A*F, G)
+
+    assert is_linear_expression(f(F,G), (F,))
+    assert is_linear_expression(f(F,G), (G,))
+    # ...
+
+    # ...
+    f = lambda F,G: cross(F, A*G)
+
+    assert is_linear_expression(f(F,G), (F,))
+    assert is_linear_expression(f(F,G), (G,))
+    # ...
+
+    # ...
+    f = lambda F,G: cross(A*F, A*G)
+
+    assert is_linear_expression(f(F,G), (F,))
+    assert is_linear_expression(f(F,G), (G,))
+    # ...
+
+    # ...
+    I = SympdeMatrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]], name='I')
+    epsilon = lambda w: 0.5*(grad(w) + Transpose(grad(w)))
+    sigma   = lambda w: kappa * div(w) * I + 2 * mu * epsilon(w)
+    f       = lambda u,v: inner(sigma(u), epsilon(v))
+
+    assert is_linear_expression(f(F,G), (F,))
+    assert is_linear_expression(f(F,G), (G,))
+    # ...
+
 #==============================================================================
 # CLEAN UP SYMPY NAMESPACE
 #==============================================================================
diff --git a/sympde/expr/tests/test_tensor_2d.py b/sympde/expr/tests/test_tensor_2d.py
index 0fc9e31d..33764286 100644
--- a/sympde/expr/tests/test_tensor_2d.py
+++ b/sympde/expr/tests/test_tensor_2d.py
@@ -2,9 +2,9 @@
 
 # TODO: - add assert to every test
 
-from sympy.core.containers import Tuple
 
 from sympde.core     import Constant
+from sympde.core     import Vector
 from sympde.calculus import grad, dot, curl, div
 #from sympde.calculus import laplace
 #from sympde.topology import dx
@@ -112,7 +112,7 @@ def test_tensorize_2d_3():
 
     bx = Constant('bx')
     by = Constant('by')
-    b  = Tuple(bx, by)
+    b = Vector([bx,by], name='b')
 
     expr = integral(domain, dot(b, grad(v)) * dot(b, grad(u)))
     a = BilinearForm((u,v), expr)
diff --git a/sympde/exterior/tests/test_datatype.py b/sympde/exterior/tests/test_datatype.py
index 9bea28d1..6a904464 100644
--- a/sympde/exterior/tests/test_datatype.py
+++ b/sympde/exterior/tests/test_datatype.py
@@ -2,7 +2,6 @@
 
 
 from sympy import symbols
-from sympy import Tuple
 from sympy import Matrix
 from sympy import srepr
 from sympy import Symbol
diff --git a/sympde/exterior/tests/test_exterior.py b/sympde/exterior/tests/test_exterior.py
index 07d28c88..caa882df 100644
--- a/sympde/exterior/tests/test_exterior.py
+++ b/sympde/exterior/tests/test_exterior.py
@@ -2,7 +2,6 @@
 
 
 from sympy import symbols
-from sympy import Tuple
 from sympy import Matrix
 from sympy import srepr
 from sympy import Symbol
diff --git a/sympde/exterior/tests/test_inference.py b/sympde/exterior/tests/test_inference.py
index 850cdaaa..ca12b2e3 100644
--- a/sympde/exterior/tests/test_inference.py
+++ b/sympde/exterior/tests/test_inference.py
@@ -2,7 +2,6 @@
 
 
 from sympy import symbols
-from sympy import Tuple
 from sympy import Matrix
 from sympy import srepr
 from sympy import Symbol
diff --git a/sympde/topology/derivatives.py b/sympde/topology/derivatives.py
index f5d3537a..14169e7d 100644
--- a/sympde/topology/derivatives.py
+++ b/sympde/topology/derivatives.py
@@ -23,6 +23,8 @@
 from sympde.core.basic    import CalculusFunction
 from sympde.core.basic    import _coeffs_registery
 from sympde.core.basic    import BasicMapping
+from sympde.core.matrices import Vector
+from sympde.core.matrices import Matrix as SympdeMatrix
 from sympde.core.algebra  import LinearOperator
 from sympde.calculus.core import minus, plus
 from sympde.calculus.core import has
@@ -82,6 +84,10 @@ def eval(cls, expr):
             args = [cls(i, evaluate=True) for i in expr]
             args = Tuple(*args)
             return Matrix([args])
+        elif isinstance(expr, Vector):
+            args = [[cls(a, evaluate=True)] for a in expr.args[0]]
+            return ImmutableDenseMatrix(args)
+
         elif isinstance(expr, (Matrix, ImmutableDenseMatrix)):
             newexpr = Matrix.zeros(*expr.shape)
             for i in range(expr.shape[0]):
diff --git a/sympde/topology/mapping.py b/sympde/topology/mapping.py
index 5c8b6589..75ad3867 100644
--- a/sympde/topology/mapping.py
+++ b/sympde/topology/mapping.py
@@ -19,12 +19,14 @@
 from sympde.core.basic        import BasicMapping
 from sympde.core.basic        import CalculusFunction
 from sympde.core.basic        import _coeffs_registery
+from sympde.core.matrices     import MatrixSymbolicExpr, MatrixElement, SymbolicTrace, Inverse
+from sympde.core.matrices     import SymbolicDeterminant, Transpose
+from sympde.core.matrices     import Vector
+from sympde.core.matrices     import Matrix as SympdeMatrix
 from sympde.calculus.core     import PlusInterfaceOperator, MinusInterfaceOperator
 from sympde.calculus.core     import grad, div, curl, laplace #, hessian
 from sympde.calculus.core     import dot, inner, outer, _diff_ops
 from sympde.calculus.core     import has, DiffOperator
-from sympde.calculus.matrices import MatrixSymbolicExpr, MatrixElement, SymbolicTrace, Inverse
-from sympde.calculus.matrices import SymbolicDeterminant, Transpose
 
 from .basic       import BasicDomain, Union, InteriorDomain
 from .basic       import Boundary, Connectivity, Interface
@@ -797,7 +799,10 @@ def eval(cls, F, v):
             the covariant transformation
         """
 
-        if not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)):
+        if isinstance(v, (Vector, SympdeMatrix)):
+            v = v.to_sympy()
+
+        elif not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)):
             raise TypeError('> Expecting a tuple, list, Tuple, Matrix')
 
         assert F.pdim == F.ldim
@@ -849,7 +854,10 @@ def eval(cls, F, v):
         if not isinstance(F, Mapping):
             raise TypeError('> Expecting a Mapping')
 
-        if not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)):
+        if isinstance(v, (Vector, SympdeMatrix)):
+            v = v.to_sympy()
+
+        elif not isinstance(v, (tuple, list, Tuple, ImmutableDenseMatrix, Matrix)):
             raise TypeError('> Expecting a tuple, list, Tuple, Matrix')
 
         M = Jacobian(F)
@@ -976,7 +984,7 @@ def eval(cls, expr, domain, **options):
         elif isinstance(expr, Transpose):
             arg = cls(expr.arg, domain)
             return Transpose(arg)
-            
+
         elif isinstance(expr, grad):
             arg = expr.args[0]
             if isinstance(mapping, InterfaceMapping):
@@ -1117,7 +1125,7 @@ def eval(cls, expr, domain, **options):
 
             elif dim == 3:
                 lgrad_arg = LogicalGrad_3d(arg)
-            
+
             grad_arg = Covariant(mapping, lgrad_arg)
             expr = grad_arg[0]
             return expr
@@ -1243,7 +1251,7 @@ def eval(cls, expr, domain, **options):
             domain  = domain.logical_domain
             det     = TerminalExpr(sqrt((J.T*J).det()), domain=domain)
             return DomainExpression(domain, ImmutableDenseMatrix([[newexpr*det]]))
-            
+
         elif isinstance(expr, Function):
             args = [cls.eval(a, domain) for a in expr.args]
             return type(expr)(*args)
diff --git a/sympde/topology/tests/test_derivatives.py b/sympde/topology/tests/test_derivatives.py
index 3694adfe..b7cf31d0 100644
--- a/sympde/topology/tests/test_derivatives.py
+++ b/sympde/topology/tests/test_derivatives.py
@@ -1,11 +1,11 @@
 # coding: utf-8
 
 from sympy import symbols
-from sympy import Tuple
 from sympy import Matrix
 from sympy import srepr
 
 from sympde.core import Constant
+from sympde.core import Vector
 from sympde.calculus import grad, dot, inner
 from sympde.topology import Domain, element_of
 from sympde.topology import get_index_derivatives_atom
@@ -25,6 +25,7 @@ def indices_as_str(a):
 
 
 # ...
+# TODO ARA add more tests for sympde/Vector and Matrix
 def test_partial_derivatives_1():
     print('============ test_partial_derivatives_1 ==============')
 
@@ -39,7 +40,7 @@ def test_partial_derivatives_1():
     V = ScalarFunctionSpace('V', mapped_domain)
 
     F,u,v,w = [element_of(V, name=i) for i in ['F', 'u', 'v', 'w']]
-    uvw = Tuple(u,v,w)
+    uvw = Vector([u,v,w], name='uvw')
 
     alpha = Constant('alpha')
     beta = Constant('beta')
@@ -54,14 +55,14 @@ def test_partial_derivatives_1():
     assert(dz(y**2) == 0)
 
     assert(dx(x*F) == F + x*dx(F))
-    assert(dx(uvw) == Matrix([[dx(u), dx(v), dx(w)]]))
-    assert(dx(uvw) + dy(uvw) == Matrix([[dx(u) + dy(u),
-                                         dx(v) + dy(v),
-                                         dx(w) + dy(w)]]))
-
-    expected = Matrix([[alpha*dx(u) + beta*dy(u),
-                        alpha*dx(v) + beta*dy(v),
-                        alpha*dx(w) + beta*dy(w)]])
+    assert(dx(uvw) == Matrix([[dx(u)], [dx(v)], [dx(w)]]))
+    assert(dx(uvw) + dy(uvw) == Matrix([[dx(u) + dy(u)],
+                                        [dx(v) + dy(v)],
+                                        [dx(w) + dy(w)]]))
+
+    expected = Matrix([[alpha*dx(u) + beta*dy(u)],
+                       [alpha*dx(v) + beta*dy(v)],
+                       [alpha*dx(w) + beta*dy(w)]])
     assert(alpha * dx(uvw) + beta * dy(uvw) == expected)
     # ...
 
diff --git a/sympde/topology/tests/test_gallery.py b/sympde/topology/tests/test_gallery.py
index 6222c9ce..bf2875c1 100644
--- a/sympde/topology/tests/test_gallery.py
+++ b/sympde/topology/tests/test_gallery.py
@@ -1,7 +1,6 @@
 # coding: utf-8
 
 from sympy.abc import x,y,z
-from sympy import Tuple
 from sympy import symbols
 
 x1, x2, x3 = symbols('x1, x2, x3')
diff --git a/sympde/topology/tests/test_mapping.py b/sympde/topology/tests/test_mapping.py
index 84ede14e..30e09df5 100644
--- a/sympde/topology/tests/test_mapping.py
+++ b/sympde/topology/tests/test_mapping.py
@@ -5,12 +5,13 @@
 from sympy.tensor import IndexedBase
 from sympy import symbols, simplify
 
+from sympde.core.matrices import Vector
 from sympde.topology import Mapping, MappedDomain
 from sympde.topology import dx, dy, dz
 from sympde.topology import dx1, dx2, dx3
 from sympde.topology import Domain
-
 from sympde.topology.mapping import Jacobian, Covariant, Contravariant
+
 # ...
 def test_mapping_1d():
     print('============ test_mapping_1d ==============')
@@ -41,7 +42,7 @@ def test_mapping_2d():
     F = Mapping('F', dim=dim)
 
     a,b = symbols('a b')
-    ab = Tuple(a, b)
+    ab = Vector([a, b], name='ab')
 
     assert(F.name == 'F')
 
@@ -81,7 +82,7 @@ def test_mapping_3d():
     F = Mapping('F', dim=dim)
 
     a,b,c = symbols('a b c')
-    abc = Tuple(a, b, c)
+    abc = Vector([a, b, c], name='abc')
 
     assert(F.name == 'F')