def METHOD_NAME():
    """
    This example illustrates how to set up a diffusion/transport model with
    a simple reaction-diffusion system in a tapering cylinder:
    | Molecule **a** diffuses with diffConst of 10e-12 m^2/s.
    | Molecule **b** diffuses with diffConst of 5e-12 m^2/s.
    | Molecule **b** also undergoes motor transport with a rate of 10e-6 m/s
    |   Thus it 'piles up' at the end of the cylinder.
    | Molecule **c** does not move: diffConst = 0.0
    | Molecule **d** does not move: diffConst = 10.0e-12 but it is buffered.
    |   Because it is buffered, it is treated as non-diffusing.
    All molecules other than **d** start out only in the leftmost (first)
    voxel, with a concentration of 1 mM. **d** is present throughout
    at 0.2 mM, except in the last voxel, where it is at 1.0 mM.
    The cylinder has a starting radius of 2 microns, and end radius of
    1 micron. So when the molecule undergoing motor transport gets to the
    narrower end, its concentration goes up.
    There is a little reaction in all compartments: ``b + d <===> c``
    As there is a high concentration of **d** in the last compartment,
    when the molecule **b** reaches the end of the cylinder, the reaction
    produces lots of **c**.
    Note that molecule **a** does not participate in this reaction.
    The concentrations of all molecules are displayed in an animation.
    """
    runtime = 20.0
    diffdt = 0.005
    plotdt = 0.1
    makeModel()
    # Set up clocks. The dsolver to know before assigning stoich
    moose.setClock( 10, diffdt ) # 10 is the standard clock for Dsolve.
    moose.setClock( 16, plotdt ) # 16 is the standard clock for Ksolve.
    a = moose.element( '/model/compartment/a' )
    b = moose.element( '/model/compartment/b' )
    c = moose.element( '/model/compartment/c' )
    d = moose.element( '/model/compartment/d' )
    moose.reinit()
    atot = sum( a.vec.n )
    btot = sum( b.vec.n )
    ctot = sum( c.vec.n )
    dtot = sum( d.vec.n )
    for t in np.arange( 0, runtime, plotdt ):
        moose.start( plotdt )
    atot2 = sum( a.vec.n )
    btot2 = sum( b.vec.n )
    ctot2 = sum( c.vec.n )
    dtot2 = sum( d.vec.n )
    msg = 'Ratio of initial to final total numbers of '
    got = np.array((atot2/atot, btot2/btot, ctot2/ctot, dtot2/dtot))
    msg += 'a=%f b=%f, c=%f, d=%f' % (tuple(got))
    print(msg)
    print('Initial to final (b+c)=%f' % (float(btot2 + ctot2) / (btot + ctot )))
    expected = np.array((1.00003087, 1.39036644, 0.92191184, 1.11427514))
    error = got - expected
    rerror = np.abs( error ) / expected
    assert np.allclose(got, expected, atol=1e-3), "Got %s, expected %s" % (got, expected)