title: meshgrid doesn't work when function contains matrices
body: <p>I wanted to plot a surface of a real function of two variables, which
  contain matrices, but otherwise the function is real. When I try to run
  the code:</p>
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  from scipy.linalg import expm, norm
  from mpl_toolkits.mplot3d import Axes3D
  from matplotlib import cm
  
  def f(x,y):
    A = np.array([[0,-1,0],[1,0,0],[0,0,0]])
    v = np.array([[1,2,3]])
    return norm(x*y*A@v.T)
  
  x = np.linspace(-10,10)
  xx, yy = np.meshgrid(x,x)
  ff = f(xx, yy)
  
  plt.figure()
  ax = fig.gca(projection='3d')
  surf = ax.plot_surface(xx, yy, ff, cmap=cm.viridis)
  </code></pre>
  <p>I always get:</p>
  <pre><code>[...]
  ValueError: operands could not be broadcast together with shapes (50,50) (3,3)
  </code></pre>
  <p>Despite all operands are taken care with respect to shape.
  What do I have to modify in order that it runs?
  Many thanks in advance</p>