title: Why is my 2D interpolant generating a matrix with swapped axes in SciPy?
body: <p>I solve a differential equation with vector inputs</p>
  
  <p><strong>y' = f(t,y), y(t_0) = y_0</strong></p>
  
  <p>where <strong>y0 = y(x)</strong></p>
  
  <p>using the explicit Euler method, which says that</p>
  
  <p><strong>y_(i+1) = y_i + h*f(t_i, y_i)</strong></p>
  
  <p>where <strong>t</strong> is a time vector, <strong>h</strong> is the step size, and <strong>f</strong> is the right-hand side of the differential equation.</p>
  
  <p>The python code for the method looks like this:</p>
  
  <pre><code>for i in np.arange(0,n-1):
          y[i+1,...] = y[i,...] + dt*myode(t[i],y[i,...])
  </code></pre>
  
  <p>The result is a <strong>k,m</strong> matrix <strong>y</strong>, where <strong>k</strong> is the size of the <strong>t</strong> dimension, and <strong>m</strong> is the size of <strong>y</strong>.</p>
  
  <p>The vectors <strong>y</strong> and <strong>t</strong> are returned.</p>
  
  <p><strong>t</strong>, <strong>x</strong>, and <strong>y</strong> are passed to <code>scipy.interpolate.RectBivariateSpline(t, x, y, kx=1, ky=1)</code>:</p>
  
  <pre><code>g = scipy.interpolate.RectBivariateSpline(t, x, y, kx=1, ky=1)
  </code></pre>
  
  <p>The resulting object <code>g</code> takes new vectors <code>ti,xi ( g(p,q) )</code> to give <code>y_int</code>, which is <strong>y</strong> interpolated at the points defined by <strong>ti</strong> and <strong>xi</strong>.</p>
  
  <p>Here is my problem:</p>
  
  <ul>
  <li><p>The documentation for RectBivariateSpline describes the <code>__call__</code> method in terms of x and y:
  <code>__call__(x, y[, mth])</code> <em>Evaluate spline at the grid points defined by the coordinate arrays</em></p></li>
  <li><p>The matplotlib documentation for plot_surface uses similar notation:
  <code>Axes3D.plot_surface(X, Y, Z, *args, **kwargs)</code>
  with the important difference that X and Y are 2D arrays which are generated by <code>numpy.meshgrid()</code>.</p></li>
  </ul>
  
  <p>When I compute simple examples, the input order is the same in both and the result is exactly what I would expect. In my explicit Euler example, however, the initial order is <code>ti,xi</code>, yet the surface plot of the interpolant output only makes sense if I reverse the order of the inputs, like so:</p>
  
  <pre><code>ax2.plot_surface(xi, ti, u, cmap=cm.coolwarm)
  </code></pre>
  
  <p>While I am glad that it works, I'm not satisfied because <strong>I cannot explain why</strong>, nor why (apart from the array geometry) it is necessary to swap the inputs. Ideally, I would like to restructure the code so that the input order is consistent.</p>
  
  <p>Here is a working code example to illustrate what I mean:</p>
  
  <pre><code># Heat equation example with explicit Euler method
  import numpy as np
  import matplotlib.pyplot as mplot
  import matplotlib.cm as cm
  import scipy.sparse as sp
  import scipy.interpolate as interp
  from mpl_toolkits.mplot3d import Axes3D
  import pdb
  
  # explicit Euler method
  def eev(myode,tspan,y0,dt):
      # Preprocessing
      # Time steps
      tspan[1] = tspan[1] + dt
      t = np.arange(tspan[0],tspan[1],dt,dtype=float)
      n = t.size
      m = y0.shape[0]
      y = np.zeros((n,m),dtype=float)
      y[0,:] = y0
  
      # explicit Euler recurrence relation
      for i in np.arange(0,n-1):
          y[i+1,...] = y[i,...] + dt*myode(t[i],y[i,...])
  
      return y,t
  
  # generate matrix A
  # u'(t) = A*u(t) + g*u(t)
  def a_matrix(n):
      aa = sp.diags([1, -2, 1],[-1,0,1],(n,n))
      return aa
  
  # System of ODEs with finite differences
  def f(t,u):
      dydt = np.divide(1,h**2)*A.dot(u)
      return dydt
  
  # homogenous Dirichlet boundary conditions
  def rbd(t):
      ul = np.zeros((t,1))
      return ul
  
  # Initial value problem -----------
  
  def main():
      # Metal rod 
      # spatial discretization
      # number of inner nodes
      m = 20
      x0 = 0
      xn = 1
      x = np.linspace(x0,xn,m+2)
      # Step size
      global h
      h = x[1]-x[0]
  
      # Initial values
      u0 = np.sin(np.pi*x)
  
      # A matrix
      global A
      A = a_matrix(m)
  
      # Time
      t0 = 0
      tend = 0.2
      # Time step width
      dt = 0.0001
      tspan = [t0,tend]
  
      # Test r for stability
      r = np.divide(dt,h**2)
      if r &lt;= 0.5:
          u,t = eev(f,tspan,u0[1:-1],dt)      
      else:
          print('r = ',r)
          print('r &gt; 0.5. Explicit Euler method will not be stable.')
  
      # Add boundary values back
      rb = rbd(t.size)
      u = np.hstack((rb,u,rb))
  
      # Interpolate heat values
      # Create interpolant. Note the parameter order
      fi = interp.RectBivariateSpline(t, x, u, kx=1, ky=1)
  
      # Create vectors for interpolant
      xi = np.linspace(x[0],x[-1],100)
      ti = np.linspace(t0,tend,100)
  
      # Compute function values from interpolant
      u_int = fi(ti,xi)
  
      # Change xi, ti in to 2D arrays
      xi,ti = np.meshgrid(xi,ti)
  
      # Create figure and axes objects
      fig3 = mplot.figure(1)
      ax3 = fig3.gca(projection='3d')
      print('xi.shape =',xi.shape,'ti.shape =',ti.shape,'u_int.shape =',u_int.shape)
  
      # Plot surface. Note the parameter order, compare with interpolant!
      ax3.plot_surface(xi, ti, u_int, cmap=cm.coolwarm)
      ax3.set_xlabel('xi')
      ax3.set_ylabel('ti')
  
  main()
  mplot.show()
  </code></pre>