title: Sympy notation for matrix multiplication
body: <p>In numpy it is standard to define matrices and dot product as shown below</p>
  
  <pre><code>import numpy as np
  a = np.array([[1, 2, 3], [4, 5, 6]])
  b = np.array([[7, 8, 9, 10], [11, 12, 13, 14], [15, 16, 17, 18]])
  print(a.shape)
  print(b.shape)
  print(a.dot(b).shape)
  </code></pre>
  
  <p>which outputs as expected:</p>
  
  <pre><code>(2, 3)
  (3, 4)
  (2, 4)
  </code></pre>
  
  <p>But to my surprise the following fails in sympy</p>
  
  <pre><code>import sympy as sp
  a = sp.Matrix([[1, 2, 3], [4, 5, 6]])
  b = sp.Matrix([[7, 8, 9, 10], [11, 12, 13, 14], [15, 16, 17, 18]])
  print(a.shape)
  print(b.shape)
  print(a.dot(b).shape)
  </code></pre>
  
  <p>which outputs</p>
  
  <pre><code>(2, 3)
  (3, 4)
  
  ---------------------------------------------------------------------------
  ShapeError                                Traceback (most recent call last)
  &lt;ipython-input-30-50c934c7fbaf&gt; in &lt;module&gt;()
        4 print(a.shape)
        5 print(b.shape)
  ----&gt; 6 print(a.dot(b).shape)
  
  ~/miniconda3/lib/python3.6/site-packages/sympy/matrices/matrices.py in dot(self, b)
     2389                 mat = mat.T
     2390                 b = b.T
  -&gt; 2391             prod = flatten((mat * b).tolist())
     2392             if len(prod) == 1:
     2393                 return prod[0]
  
  ~/miniconda3/lib/python3.6/site-packages/sympy/core/decorators.py in binary_op_wrapper(self, other)
      130                     else:
      131                         return f(self)
  --&gt; 132             return func(self, other)
      133         return binary_op_wrapper
      134     return priority_decorator
  
  ~/miniconda3/lib/python3.6/site-packages/sympy/matrices/common.py in __mul__(self, other)
     2006             if self.shape[1] != other.shape[0]:
     2007                 raise ShapeError("Matrix size mismatch: %s * %s." % (
  -&gt; 2008                     self.shape, other.shape))
     2009 
     2010         # honest sympy matrices defer to their class's routine
  
  ShapeError: Matrix size mismatch: (3, 2) * (4, 3).
  </code></pre>
  
  <p>This is confusing for me!</p>
  
  <p>Why the inconsistency between numpy and sympy?</p>
  
  <p>Why can't I find a warning of this behaviour in sympy's documentation?</p>
  
  <p>How can I correctly compute the dot product of two matrices in sympy?</p>
  
  <p>May I suggest that the documentation of sympy has an entry on the differences of syntax between numpy and sympy. (I would gladly contribute but I have no idea on those differences)</p>