Python vector field of ODE of three variables
Question: I am trying to plot a vector field of a ODE model with three variables. I
would like to average the vectors along the third axis, and present the vector
field together with the information of the standard deviation of their values.
The ODE system is:

    
    
    a  = 1.
    b1 = 0.1
    b2 = 0.11
    c1 = 1.5
    c2 = 1.6
    d  = 0.75
    def dudt(a,b1,b2,u,v1,v2):
        return a*u - b1*u*v1 - b2*u*v2
    def dv1dt(d,c1,b1,u,v1):
        return -c1*v1 + d*b1*u*v1
    def dv2dt(d,c2,b2,u,v2):
        return -c2*v2 + d*b2*u*v2
    

The function that I am currently using is:

    
    
    import numpy as np
    import matplotlib.pyplot as plt
    def plotVF(mS=None, density= 20,color='k'):
        mB1 = np.linspace(0,1.1,int(density))
        mB2 = np.linspace(0,1.1,int(density))
        if mS==None:
            mS  = np.linspace(0,1.1,int(density))
            B1,B2,S = np.meshgrid(mB1,mB2,mS)
            average=True
        else:
            B1,B2 = np.meshgrid(mB1,mB2)
            S = mS
            average=False
        DB1 = dv1dt(d,c1,b1,S,B1)
        DB2 = dv2dt(d,c2,b2,S,B2)
        DS  = dudt(a,b1,b2,S,B1,B2)
        if average:
            print "Averaging"
            DB1std = np.std(DB1,axis=2)
            DB2std = np.std(DB2,axis=2)
            DB1 = np.mean(DB1,axis=2)
            DB2 = np.mean(DB2,axis=2)
            DS  = np.mean(DS,axis=2)
            vecstd = np.hypot(DB1std,DB2std)
            plt.imshow(vecstd)
            plt.colorbar()
            B1,B2 = np.meshgrid(mB1,mB2)
        M = (np.hypot(DB1, DB2, DS))
        M[ M == 0] = 1.
        DB1=DB1/M
        DB2=DB2/M
        DS=DS/M
        print B1.shape,B2.shape,DB1.shape,DB2.shape
        plt.quiver(B1, B2, DB1, DB2, pivot='mid', color=color)
        plt.xlim(0,1.1), plt.ylim(0,1.1)
        plt.grid('on')
        plt.show()
    

It gives me that the standard deviation along the third axis is zero, which
does not make sense.
[![Vector_Field](http://i.stack.imgur.com/0HSyY.png)](http://i.stack.imgur.com/0HSyY.png)
Someone has an idea what am I doing wrong?


Answer: Your code is almost perfectly fine. There's only one problem: you're plotting
the colormap with a vanilla call to `plt.imshow`.

As its name suggests, `imshow` is used for plotting images. As such, it by
default doesn't expect coordinate inputs, just a single array containing the
pixel data. This implies that a simple call to `imshow` will have axes limits
corresponding to the number of pixels in your image -- in your case the
dimensions of your 2d data arrays. If you take a look at the image directly
created by `imshow`, you'll see that the limits go up to `x,y=20`. Later you
set new limits according to your actual underlying mesh, truncating your plot
to the first 2 data points.

The solution is to explicitly tell `plt.imshow()` where you want your plot to
reside in coordinate space:

    
    
    plt.imshow(vecstd, extent=[B1.min(),B1.max(),B2.min(),B2.max()], origin='lower')
    

The first keyword argument `extent` gives the `x` and `y` limits into which
the data should be plotted. Note the important second keyword argument,
`origin`. By default `imshow` plots things "upside down" in order to _not_
plot actual images upside down. When you're using `imshow` to plot stuff
defined with Cartesian coordinates, you have to tell it that the origin of the
coordinate system should not be the _upper_ left corner of the figure (as for
images), bur rather the _lower_ left corner (as for regular plots).

