Starting with the advection equation is conservative form,

$$
\frac{\partial u}{\partial t} = -\frac{\partial (\boldsymbol{v} u)}{\partial x} + s(x,t)
$$

The Crank-Nicolson method consists of a time averaged centered difference.

$$\frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} = -\boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( u_{j+1}^{n} - u_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( u_{j+1}^{n+1} - u_{j-1}^{n+1} \right) \right] + s(x,t)$$

Regarding the notation, *subscripts* are for points in space, and the *superscripts* are for points in time.

The points at $n+1$ are in the future: they are unknowns. We now need to rearrange the above equation so that all *knowns* are on the r.h.s and *unknowns* are on the l.h.s..

Making the substitution,

$$
r = \frac{\boldsymbol{v}}{2}\frac{\Delta t}{\Delta x}
$$

gives,

$$-\beta r\phi_{j-1}^{n+1} + \phi_{j}^{n+1} + \beta r\phi_{j+1}^{n+1} =  (1-\beta)r\phi_{j-1}^{n} + \phi_{j}^{n} - (1-\beta)r\phi_{j+1}^{n}$$

This is the advection equation discretised using the Crank-Nicolson method. You can write it as a matrix equation,

$$
    \begin{pmatrix}
      1                 & \beta r        	&     				&       	&          0		\\
      -\beta r    		& 1         		& \beta r       	&     	    &		\\
              			&  \ddots   		&   \ddots  		& \ddots	&            		\\
						&                   & -\beta r  		    & 1         & \beta r 	        \\
                 0      &                   &      				& -\beta r  	& 1          		\\
   \end{pmatrix}
   \begin{pmatrix}
u_1^{n+1} \\
u_2^{n+1} \\
\vdots \\
u_{J-1}^{n+1} \\
u_{J}^{n+1} \\
   \end{pmatrix} = 
    \begin{pmatrix}
      1                    & -(1 - \beta)r        		&     						&       	          &   0		 	        \\
      (1 - \beta)r        & 1 			                & -(1 - \beta)r 		        &		     	      &                     \\
              			   &  \ddots   					&   \ddots  				& \ddots			  &                     \\
						   &                            & (1 - \beta)r  		    & 1	                  & -(1 - \beta)r	\\
                 0         &                            &      						&(1 - \beta)r        & 1	             \\
   \end{pmatrix}   \begin{pmatrix}
u_1^{n} \\
u_2^{n} \\
\vdots \\
u_{J-1}^{n} \\
u_{J}^{n} \\
   \end{pmatrix} 
$$
Setting $\beta=1/2$ will give you trapezoidal integration in time, so for Crank-Nicolson this is what you want.

A few words of warning. This is basic solution you wanted, but you will need to include some sort of **boundary condition** for a well-posed problem. Also,
Crank-Nicolson is not necessarily the best method for the advection equation. It is second order accurate and [unconditionally stable][1], which is fantastic. However it will generate (as with all centered difference stencils) spurious oscillation if you have very sharp peaked solutions or initial conditions.

I wrote the following code for you in Python, it should get you started. The code solves the advection equation for an initial Gaussian curve moving to the right with constant velocity.

![Gaussian curve moving to the right with constant velocity][2]

    from __future__ import division
    from scipy.sparse import spdiags
    from scipy.sparse.linalg import spsolve
    import numpy as np
    import pylab
    
    def make_advection_matrices(z, r):
        """Return matrices A and M for advection equations"""
        ones = np.ones(len(z))
        A = spdiags( [-beta*r, ones, beta*r], (-1,0,1), len(z), len(z) )
        M = spdiags( [(1-beta) * r, ones, -(1-beta) * r], (-1,0,1), len(z), len(z) )
        return A.tocsr(), M.tocsr()
    
    def plot_iteration(z, u, iteration):
        """Plot the solver progress"""
        pylab.plot(z, u, label="Iteration %d" % iteration)
    
    # Set up basic constants
    beta = 0.5
    J = 200 # total number of mesh points
    z = np.linspace(-10,10,J) # vertices
    dz = abs(z[1]-z[0]) # space step
    dt = 0.2    # time step
    v = 2 * np.ones(len(z)) # velocity field (constant)
    r = v / 2 * dt / dz
    
    # Initial conditions (peak function)
    gaussian = lambda z, height, position, hwhm: height * np.exp(-np.log(2) * ((z - position)/hwhm)**2)
    u_init = gaussian(z, 1, -3, 2)
    
    A, M = make_advection_matrices(z, r)
    u = u_init
    for i in range(10):
        u = spsolve(A, M * u)
        plot_iteration(z, u, i)
    
    pylab.legend()
    pylab.show()


  [1]: http://scicomp.stackexchange.com/questions/5402/is-crank-nicolson-a-stable-discretization-scheme-for-reaction-diffusion-advectio
  [2]: https://i.sstatic.net/zFaVE.gif