diff --git a/demos/feature_demo.py b/demos/feature_demo.py
new file mode 100644
index 00000000..78a86995
--- /dev/null
+++ b/demos/feature_demo.py
@@ -0,0 +1,126 @@
+# This script creates a set of PDFs that illustrate the effect on the SED of
+# successively turning on various options or changing the value of some
+# variables.
+
+import numpy as np
+import matplotlib.pyplot as pl
+from matplotlib.backends.backend_pdf import PdfPages
+
+import fsps
+
+
+def makefig(sps, tage=13.7, oldspec=None, **plotkwargs):
+    w, spec = sps.get_spectrum(tage=tage)
+    fig, ax = pl.subplots()
+    if oldspec is not None:
+        ax.plot(w, oldspec / w * 1e19, color='gray', linewidth=2, alpha=0.5)
+    ax.plot(w, spec / w * 1e19, 'C2', linewidth=2)
+    return fig, ax, spec
+
+
+def prettify(fig, ax, label=None):
+    ax.set_xlim(0.9e3, 1e6)
+    ax.set_xscale('log')
+    ax.set_ylim(0.01, 2)
+    #ax.set_yscale('log')
+    ax.set_xlabel('rest-frame $\lambda$ ($\AA$)', fontsize=20)
+    ax.set_ylabel('$\lambda \, f_\lambda$', fontsize=20)
+    ax.tick_params(axis='both', which='major', labelsize=16)
+    if label is not None:
+        ax.text(0.63, 0.85, label, transform=ax.transAxes, fontsize=16)
+    
+    return fig, ax
+
+if __name__ == "__main__":
+
+    pl.rc('text', usetex=True)
+    pl.rc('font', family='serif')
+    pl.rc('axes', grid=False)
+    pl.rc('xtick', direction='in')
+    pl.rc('ytick', direction='in')
+    pl.rc('xtick', top=True)
+    pl.rc('ytick', right=True)
+
+    sps = fsps.StellarPopulation(zcontinuous=1)
+    pdf = PdfPages('features.pdf')
+
+    # Basic spectrum
+    sps.params['sfh'] = 4
+    sps.params['tau'] = 5.0
+    sps.params['logzsol'] = 0.0
+    sps.params['dust_type'] = 4 # kriek and Conroy
+    sps.params['imf_type'] = 2  #kroupa
+    sps.params['imf3'] = 2.3
+    fig, ax, spec = makefig(sps)
+    fig, ax = prettify(fig, ax, label="$\\tau=5$, Age$=13.7$,\n$\log Z/Z_\odot=0.0$")
+    pdf.savefig(fig)
+    pl.close(fig)
+    
+    # change IMF
+    sps.params['imf3'] = 2.5
+    fig, ax, spec = makefig(sps, oldspec=spec)
+    fig, ax = prettify(fig, ax, label="IMF slope")
+    pdf.savefig(fig)
+    
+    # Attenuate
+    sps.params['add_dust_emission'] = False
+    sps.params['dust2'] = 0.2
+    fig, ax, spec = makefig(sps, oldspec=spec)
+    fig, ax = prettify(fig, ax, label="Dust Attenuation")
+    pdf.savefig(fig)
+    pl.close(fig)
+    
+    # Dust emission
+    sps.params['add_dust_emission'] = True
+    fig, ax, spec = makefig(sps, oldspec=spec)
+    fig, ax = prettify(fig, ax, label="Dust Emission")
+    pdf.savefig(fig)
+    pl.close(fig)
+
+    # Dust temperature
+    sps.params['duste_umin'] = 10
+    fig, ax, spec = makefig(sps, oldspec=spec)
+    fig, ax = prettify(fig, ax, label="Dust SED\n(Draine)")
+    pdf.savefig(fig)
+    pl.close(fig)
+    
+    # AGN emission
+    sps.params['fagn'] = 0.3
+    fig, ax, spec = makefig(sps, oldspec=spec)
+    fig, ax = prettify(fig, ax, label="AGN dust\n(Nenkova)")
+    pdf.savefig(fig)
+    pl.close(fig)
+    
+    # Nebular emission
+    sps.params['add_neb_emission'] = True
+    sps.params['gas_logu'] = -3.5
+    fig, ax, spec = makefig(sps, oldspec=spec)
+    fig, ax = prettify(fig, ax, label="Neb. emission\n(Byler)")
+    pdf.savefig(fig)
+    pl.close(fig)
+    
+    # change logu
+    sps.params['gas_logu'] = -1.0
+    fig, ax, spec = makefig(sps, oldspec=spec)
+    fig, ax = prettify(fig, ax, label="Change U$_{neb}$")
+    pdf.savefig(fig)
+    pl.close(fig)
+    
+    # change logz
+    sps.params['logzsol'] = -0.5
+    sps.params['gas_logz'] = -0.5
+    fig, ax, spec = makefig(sps, oldspec=spec)
+    fig, ax = prettify(fig, ax, label="$\log Z/Z_\odot=-0.5$")
+    pdf.savefig(fig)
+    pl.close(fig)
+
+    # IGM absorption
+    sps.params['zred'] = 6.0
+    sps.params['add_igm_absorption'] = True
+    fig, ax, spec = makefig(sps, oldspec=spec)
+    fig, ax = prettify(fig, ax, label="IGM attenuation\n(Madau, $z=6$)")
+    pdf.savefig(fig)
+    pl.close(fig)
+    
+    pdf.close()
+    
diff --git a/docs/stellarpop_api.rst b/docs/stellarpop_api.rst
index 3b0edbb1..df3ef5ca 100644
--- a/docs/stellarpop_api.rst
+++ b/docs/stellarpop_api.rst
@@ -14,7 +14,7 @@ and some dust with a Calzetti et al (2000) extinction curve::
    >>> import fsps
    >>> sp = fsps.StellarPopulation(compute_vega_mags=False, zcontinuous=1,
                                    sfh=0, logzsol=0.0, dust_type=2, dust2=0.2)
-   >>> sp.libraries()
+   >>> sp.libraries
    ('pdva', 'miles')
 
 The last line indicates that we are using the Padova isochrones and MILES spectral library.
diff --git a/fsps/fsps.f90 b/fsps/fsps.f90
index 09690c70..888689ba 100644
--- a/fsps/fsps.f90
+++ b/fsps/fsps.f90
@@ -249,6 +249,19 @@ subroutine compute_zdep(ns,n_age,ztype)
        call compsp(0,1,outfile,mass,lbol,spec,pset,ocompsp)
     endif
 
+    if (ztype .eq. 3) then
+       ! Build the SSPs for *every* metallicity and feed all of them to compsp
+       ! for z-dependent tabular
+       do zmet=1,nz
+          if (has_ssp(zmet) .eq. 0) then
+             call ssp(zmet)
+          endif
+       enddo
+       call compsp(0,nz,outfile,mass_ssp_zz,lbol_ssp_zz,&
+            spec_ssp_zz,pset,ocompsp)
+    endif
+
+    
   end subroutine
 
   subroutine get_spec(ns,n_age,spec_out)
diff --git a/fsps/fsps.py b/fsps/fsps.py
index c253b4e1..8c87241d 100644
--- a/fsps/fsps.py
+++ b/fsps/fsps.py
@@ -50,6 +50,11 @@ class StellarPopulation(object):
         * 2: The SSPs are convolved with a metallicity distribution function
           specified by the ``logzsol`` and ``pmetals`` parameters. The value of
           ``zmet`` is ignored.
+        * 3: Use all available SSP metallicities when computing the composite
+          model, for use exclusively with tabular SFHs where the metallicity
+          evolution as function of age is given (see `set_tabular_sfh()`).  The
+          values of ``zmet`` and ``logzsol`` are ignored.  Furthermore
+          ``add_neb_emission`` must be set to False.
 
         Can only be changed during initialization.
 
@@ -837,19 +842,24 @@ def set_tabular_sfh(self, age, sfr, Z=None):
             length as ``age``.
 
         :param Z: (optional)
-            The metallicity at each age.  Currently this is ignored, and the
-            value of ``zmet`` or ``logzsol`` is used for all ages.  Thus
-            setting this parameter will result in a NotImplementedError.
+            The metallicity at each age, in units of absolute metallicity
+            (e.g. Z=0.019 for solar with the Padova isochrones and MILES
+            stellar library).
         """
         assert len(age) == len(sfr)
         ntab = len(age)
         if Z is None:
             Z = np.zeros(ntab)
+            assert self._zcontinuous != 3
         else:
-            raise(NotImplementedError)
             assert len(Z) == ntab
+            assert np.all(Z >= 0), "All values of Z must be greater than or equal 0."
+            assert self._zcontinuous == 3, "_zcontinuous must be 3 for multi-Z tabular."
+            assert self.params["add_neb_emission"] is False, ("Cannot compute nebular emission "
+                                                              "with multi-metallicity tabular SFH.")
+
         driver.set_sfh_tab(age*1e9, sfr, Z)
-        if self.params['sfh'] == 3:
+        if self.params["sfh"] == 3:
             self.params.dirtiness = max(1, self.params.dirtiness)
         else:
             print("Warning: You are setting a tabular SFH, "
diff --git a/fsps/tests.py b/fsps/tests.py
index 2a4f289d..4c300de4 100644
--- a/fsps/tests.py
+++ b/fsps/tests.py
@@ -13,6 +13,7 @@
 
 
 def _reset_default_params():
+    pop._zcontinuous = 1
     for k in pop.params.all_params:
         pop.params[k] = default_params[k]
 
@@ -145,7 +146,18 @@ def test_tabular():
     pop.params['logzsol'] = -1
     w, spec_lowz = pop.get_spectrum(tage=age.max())
     assert not np.allclose(spec / spec_lowz - 1., 0.)
-
+    pop._zcontinuous = 3
+    pop.set_tabular_sfh(age, sfr, z)
+    w, spec_multiz = pop.get_spectrum(tage=age.max())
+    assert not np.allclose(spec_lowz / spec_multiz - 1., 0.)
+    pop._zcontinuous = 1
+    pop.set_tabular_sfh(age, sfr)
+    # get mass weighted metallicity
+    mbin = np.gradient(age) * sfr
+    mwz = (z * mbin).sum() / mbin.sum()
+    pop.params['logzsol'] = np.log10(mwz/0.019)
+    w, spec_onez = pop.get_spectrum(tage=age.max())
+    assert not np.allclose(spec_onez / spec_multiz - 1., 0.)
 
 def test_mformed():
     _reset_default_params()