diff --git a/rrompy/hfengines/base/fenics_engine_base.py b/rrompy/hfengines/base/fenics_engine_base.py
index 7d46bc5..0691f33 100644
--- a/rrompy/hfengines/base/fenics_engine_base.py
+++ b/rrompy/hfengines/base/fenics_engine_base.py
@@ -1,402 +1,402 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from os import path, mkdir
 import fenics as fen
 import numpy as np
 from matplotlib import pyplot as plt
 from .numpy_engine_base import NumpyEngineBase
 from rrompy.utilities.base.types import (Np1D, strLst, FenFunc, Tuple, List,
                                          FigHandle)
 from rrompy.utilities.base.data_structures import purgeList, getNewFilename
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.solver.fenics import L2NormMatrix, fenplot, interp_project
 from .boundary_conditions import BoundaryConditions
 from rrompy.utilities.exception_manager import RROMPyException
 
 __all__ = ['FenicsEngineBase']
 
 class FenicsEngineBase(NumpyEngineBase):
     """Generic solver for parametric fenics problems."""
 
     def __init__(self, degree_threshold : int = np.inf,
                  verbosity : int = 10, timestamp : bool = True):
         super().__init__(verbosity = verbosity, timestamp = timestamp)
         self.BCManager = BoundaryConditions("Dirichlet")
         self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
         self.degree_threshold = degree_threshold
 
     @property
     def V(self):
         """Value of V."""
         return self._V
     @V.setter
     def V(self, V):
         if not type(V).__name__ == 'FunctionSpace':
             raise RROMPyException("V type not recognized.")
         self.dsToBeSet = True
         self._V = V
         self.u = fen.TrialFunction(V)
         self.v = fen.TestFunction(V)
 
     @property
     def spacedim(self):
         if hasattr(self, "_V"): return self.V.dim()
         return super().spacedim
 
     def autoSetDS(self):
         """Set FEniCS boundary measure based on boundary function handles."""
         if self.dsToBeSet:
             vbMng(self, "INIT", "Initializing boundary measures.", 20)
             mesh = self.V.mesh()
             NB = self.NeumannBoundary
             RB = self.RobinBoundary
             boundary_markers = fen.MeshFunction("size_t", mesh,
                                                 mesh.topology().dim() - 1)
             NB.mark(boundary_markers, 0)
             RB.mark(boundary_markers, 1)
             self.ds = fen.Measure("ds", domain = mesh,
                                   subdomain_data = boundary_markers)
             self.dsToBeSet = False
             vbMng(self, "DEL", "Done assembling boundary measures.", 20)
 
     def buildEnergyNormForm(self):
         """
         Build sparse matrix (in CSR format) representative of scalar product.
         """
         vbMng(self, "INIT", "Assembling energy matrix.", 20)
         self.energyNormMatrix = L2NormMatrix(self.V)
         vbMng(self, "DEL", "Done assembling energy matrix.", 20)
 
     def buildEnergyNormDualForm(self):
         """
         Build sparse matrix (in CSR format) representative of dual scalar
             product without duality.
         """
         if not hasattr(self, "energyNormMatrix"):
             self.buildEnergyNormForm()
         self.energyNormDualMatrix = self.energyNormMatrix
 
     def liftDirichletData(self) -> Np1D:
         """Lift Dirichlet datum."""
         if not hasattr(self, "_liftedDirichletDatum"):
             liftRe = interp_project(self.DirichletDatum[0], self.V)
             liftIm = interp_project(self.DirichletDatum[1], self.V)
             self._liftedDirichletDatum = (np.array(liftRe.vector())
                                   + 1.j * np.array(liftIm.vector()))
         return self._liftedDirichletDatum
 
     def reduceQuadratureDegree(self, fun:FenFunc, name:str):
         """Check whether to reduce compiler parameters to degree threshold."""
         if not np.isinf(self.degree_threshold):
             from ufl.algorithms.estimate_degrees import (
                                       estimate_total_polynomial_degree as ETPD)
             try:
                 deg = ETPD(fun)
             except:
                 return False
             if deg > self.degree_threshold:
                 vbMng(self, "MAIN", 
                       ("Reducing quadrature degree from {} to {} for "
                        "{}.").format(deg, self.degree_threshold, name), 15)
                 return True
         return False
 
     def iterReduceQuadratureDegree(self, funsNames:List[Tuple[FenFunc, str]]):
         """
         Iterate reduceQuadratureDegree over list and define reduce compiler
             parameters.
         """
         if funsNames is not None:
             for fun, name in funsNames:
                 if self.reduceQuadratureDegree(fun, name):
                     return {"quadrature_degree" : self.degree_threshold}
         return {}
 
     def plot(self, u:Np1D, warping : List[callable] = None,
              is_state : bool = False, name : str = "u", save : str = None,
              what : strLst = 'all', forceNewFile : bool = True,
              saveFormat : str = "eps", saveDPI : int = 100, show : bool = True,
              colorMap : str = "jet", fenplotArgs : dict = {},
              **figspecs) -> Tuple[FigHandle, str]:
         """
         Do some nice plots of the complex-valued function with given dofs.
 
         Args:
             u: numpy complex array with function dofs.
             warping(optional): Domain warping functions.
             is_state(optional): whether given u is value before multiplication
                 by c. Defaults to False.
             name(optional): Name to be shown as title of the plots. Defaults to
                 'u'.
             save(optional): Where to save plot(s). Defaults to None, i.e. no
                 saving.
             what(optional): Which plots to do. If list, can contain 'ABS',
                 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
                 Defaults to 'ALL'.
             forceNewFile(optional): Whether to create new output file.
             saveFormat(optional): Format for saved plot(s). Defaults to "eps".
             saveDPI(optional): DPI for saved plot(s). Defaults to 100.
             show(optional): Whether to show figure. Defaults to True.
             colorMap(optional): Pyplot colormap. Defaults to 'jet'.
             fenplotArgs(optional): Optional arguments for fenplot.
             figspecs(optional key args): Optional arguments for matplotlib
                 figure creation.
 
         Returns:
             Output filename and figure handle.
         """
         if not is_state and not self.isCEye:
             return super().plot(u, warping, False, name, save, what,
                                 forceNewFile, saveFormat, saveDPI, show,
-                                fenplotArgs, **figspecs)
+                                colorMap, fenplotArgs, **figspecs)
         if isinstance(what, (str,)):
             if what.upper() == 'ALL':
                 what = ['ABS', 'PHASE', 'REAL', 'IMAG']
             else:
                 what = [what]
         what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'],
                          listname = self.name() + ".what", baselevel = 1)
         if len(what) == 0: return
         if 'figsize' not in figspecs.keys():
             figspecs['figsize'] = plt.figaspect(1. / len(what))
 
         subplotidx = 0
         fig = plt.figure(**figspecs)
         plt.set_cmap(colorMap)
         if 'ABS' in what:
             uAb = fen.Function(self.V)
             uAb.vector().set_local(np.abs(u))
             subplotidx = subplotidx + 1
             ax = fig.add_subplot(1, len(what), subplotidx)
             p = fenplot(uAb, warping = warping, title = "|{0}|".format(name),
                         **fenplotArgs)
             if self.V.mesh().geometric_dimension() > 1:
                 fig.colorbar(p, ax = ax)
         if 'PHASE' in what:
             uPh = fen.Function(self.V)
             uPh.vector().set_local(np.angle(u))
             subplotidx = subplotidx + 1
             ax = fig.add_subplot(1, len(what), subplotidx)
             p = fenplot(uPh, warping = warping,
                         title = "phase({0})".format(name), **fenplotArgs)
             if self.V.mesh().geometric_dimension() > 1:
                 fig.colorbar(p, ax = ax)
         if 'REAL' in what:
             uRe = fen.Function(self.V)
             uRe.vector().set_local(np.real(u))
             subplotidx = subplotidx + 1
             ax = fig.add_subplot(1, len(what), subplotidx)
             p = fenplot(uRe, warping = warping, title = "Re({0})".format(name),
                         **fenplotArgs)
             if self.V.mesh().geometric_dimension() > 1:
                 fig.colorbar(p, ax = ax)
         if 'IMAG' in what:
             uIm = fen.Function(self.V)
             uIm.vector().set_local(np.imag(u))
             subplotidx = subplotidx + 1
             ax = fig.add_subplot(1, len(what), subplotidx)
             p = fenplot(uIm, warping = warping, title = "Im({0})".format(name),
                         **fenplotArgs)
             if self.V.mesh().geometric_dimension() > 1:
                 fig.colorbar(p, ax = ax)
         plt.tight_layout()
         if save is not None:
             save = save.strip()
             if forceNewFile:
                 fileOut = getNewFilename("{}_fig_".format(save), saveFormat)
             else:
                 fileOut = "{}_fig.{}".format(save, saveFormat)
             fig.savefig(fileOut, format = saveFormat, dpi = saveDPI)
         else: fileOut = None
         if show: plt.show()
         if fileOut is None: return fig
         return fig, fileOut
 
     def plotmesh(self, warping : List[callable] = None, name : str = "Mesh",
                  save : str = None, forceNewFile : bool = True,
                  saveFormat : str = "eps", saveDPI : int = 100,
                  show : bool = True, fenplotArgs : dict = {},
                  **figspecs) -> Tuple[FigHandle, str]:
         """
         Do a nice plot of the mesh.
 
         Args:
             u: numpy complex array with function dofs.
             warping(optional): Domain warping functions.
             name(optional): Name to be shown as title of the plots. Defaults to
                 'u'.
             save(optional): Where to save plot(s). Defaults to None, i.e. no
                 saving.
             forceNewFile(optional): Whether to create new output file.
             saveFormat(optional): Format for saved plot(s). Defaults to "eps".
             saveDPI(optional): DPI for saved plot(s). Defaults to 100.
             show(optional): Whether to show figure. Defaults to True.
             fenplotArgs(optional): Optional arguments for fenplot.
             figspecs(optional key args): Optional arguments for matplotlib
                 figure creation.
 
         Returns:
             Output filename and figure handle.
         """
         fig = plt.figure(**figspecs)
         fenplot(self.V.mesh(), warping = warping, **fenplotArgs)
         plt.tight_layout()
         if save is not None:
             save = save.strip()
             if forceNewFile:
                 fileOut = getNewFilename("{}_msh_".format(save), saveFormat)
             else:
                 fileOut = "{}_msh.{}".format(save, saveFormat)
             fig.savefig(fileOut, format = saveFormat, dpi = saveDPI)
         else: fileOut = None
         if show: plt.show()
         if fileOut is None: return fig
         return fig, fileOut
 
     def outParaview(self, u:Np1D, warping : List[callable] = None,
                     is_state : bool = False, name : str = "u",
                     filename : str = "out", time : float = 0.,
                     what : strLst = 'all', forceNewFile : bool = True,
                     folder : bool = False, filePW = None) -> str:
         """
         Output complex-valued function with given dofs to ParaView file.
 
         Args:
             u: numpy complex array with function dofs.
             warping(optional): Domain warping functions.
             is_state(optional): whether given u is value before multiplication
                 by c. Defaults to False.
             name(optional): Base name to be used for data output.
             filename(optional): Name of output file.
             time(optional): Timestamp.
             what(optional): Which plots to do. If list, can contain 'MESH',
                 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
                 'ALL'. Defaults to 'ALL'.
             forceNewFile(optional): Whether to create new output file.
             folder(optional): Whether to create an additional folder layer.
             filePW(optional): Fenics File entity (for time series).
 
         Returns:
             Output filename.
         """
         if not is_state and not self.isCEye:
             raise RROMPyException(("Cannot output to Paraview non-state "
                                    "object."))
         if isinstance(what, (str,)):
             if what.upper() == 'ALL':
                 what = ['MESH', 'ABS', 'PHASE', 'REAL', 'IMAG']
             else:
                 what = [what]
         what = purgeList(what, ['MESH', 'ABS', 'PHASE', 'REAL', 'IMAG'],
                          listname = self.name() + ".what", baselevel = 1)
         if len(what) == 0: return
 
         if filePW is None:
             if folder:
                 if not path.exists(filename + "/"):
                     mkdir(filename)
                 idxpath = filename.rfind("/")
                 filename += "/" + filename[idxpath + 1 :]
             if forceNewFile:
                 filePW = fen.File(getNewFilename(filename, "pvd"))
             else:
                 filePW = fen.File("{}.pvd".format(filename))
         if warping is not None:
             fen.ALE.move(self.V.mesh(),
                          interp_project(warping[0], self.V.mesh()))
         if what == ['MESH']:
             filePW << (self.V.mesh(), time)
         if 'ABS' in what:
             uAb = fen.Function(self.V, name = "{}_ABS".format(name))
             uAb.vector().set_local(np.abs(u))
             filePW << (uAb, time)
         if 'PHASE' in what:
             uPh = fen.Function(self.V, name = "{}_PHASE".format(name))
             uPh.vector().set_local(np.angle(u))
             filePW << (uPh, time)
         if 'REAL' in what:
             uRe = fen.Function(self.V, name = "{}_REAL".format(name))
             uRe.vector().set_local(np.real(u))
             filePW << (uRe, time)
         if 'IMAG' in what:
             uIm = fen.Function(self.V, name = "{}_IMAG".format(name))
             uIm.vector().set_local(np.imag(u))
             filePW << (uIm, time)
         if warping is not None:
             fen.ALE.move(self.V.mesh(),
                          interp_project(warping[1], self.V.mesh()))
         return filePW
 
     def outParaviewTimeDomain(self, u:Np1D, omega:float,
                               warping : List[callable] = None,
                               is_state : bool = False,
                               timeFinal : float = None, 
                               periodResolution : int = 20, name : str = "u",
                               filename : str = "out",
                               forceNewFile : bool = True,
                               folder : bool = False) -> str:
         """
         Output complex-valued function with given dofs to ParaView file,
             converted to time domain.
 
         Args:
             u: numpy complex array with function dofs.
             omega: frequency.
             warping(optional): Domain warping functions.
             is_state(optional): whether given u is value before multiplication
                 by c. Defaults to False.
             timeFinal(optional): final time of simulation.
             periodResolution(optional): number of time steps per period.
             name(optional): Base name to be used for data output.
             filename(optional): Name of output file.
             forceNewFile(optional): Whether to create new output file.
             folder(optional): Whether to create an additional folder layer.
 
         Returns:
             Output filename.
         """
         if not is_state and not self.isCEye:
             raise RROMPyException(("Cannot output to Paraview non-state "
                                    "object."))
         if folder:
             if not path.exists(filename + "/"):
                 mkdir(filename)
             idxpath = filename.rfind("/")
             filename += "/" + filename[idxpath + 1 :]
         if forceNewFile:
             filePW = fen.File(getNewFilename(filename, "pvd"))
         else:
             filePW = fen.File("{}.pvd".format(filename))
         omega = np.abs(omega)
         t = 0.
         dt = 2. * np.pi / omega / periodResolution
         if timeFinal is None: timeFinal = 2. * np.pi / omega - dt
         if warping is not None:
             fen.ALE.move(self.V.mesh(),
                          interp_project(warping[0], self.V.mesh()))
         for j in range(int(np.ceil(timeFinal / dt)) + 1):
             ut = fen.Function(self.V, name = name)
             ut.vector().set_local(np.real(u) * np.cos(omega * t)
                                 + np.imag(u) * np.sin(omega * t))
             filePW << (ut, t)
             t += dt
         if warping is not None:
             fen.ALE.move(self.V.mesh(),
                          interp_project(warping[1], self.V.mesh()))
         return filePW
diff --git a/rrompy/hfengines/base/numpy_engine_base.py b/rrompy/hfengines/base/numpy_engine_base.py
index 4a44511..a84f6a5 100644
--- a/rrompy/hfengines/base/numpy_engine_base.py
+++ b/rrompy/hfengines/base/numpy_engine_base.py
@@ -1,113 +1,114 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from matplotlib import pyplot as plt
 from .hfengine_base import HFEngineBase
 from rrompy.utilities.base.types import Np1D, strLst, List, Tuple, FigHandle
 from rrompy.utilities.base.data_structures import purgeList, getNewFilename
 
 __all__ = ['NumpyEngineBase']
     
 class NumpyEngineBase(HFEngineBase):
     """Generic solver for parametric matricial problems."""
 
     def plot(self, u:Np1D, warping : List[callable] = None,
              is_state : bool = False, name : str = "u", save : str = None,
              what : strLst = 'all', forceNewFile : bool = True,
              saveFormat : str = "eps", saveDPI : int = 100, show : bool = True,
              colorMap : str = "jet", pyplotArgs : dict = {},
              **figspecs) -> Tuple[FigHandle, str]:
         """
         Do some nice plots of the complex-valued function with given dofs.
 
         Args:
             u: numpy complex array with function dofs.
             name(optional): Name to be shown as title of the plots. Defaults to
                 'u'.
             is_state(optional): whether given u is value before multiplication
                 by c. Defaults to False.
             save(optional): Where to save plot(s). Defaults to None, i.e. no
                 saving.
             what(optional): Which plots to do. If list, can contain 'ABS',
                 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
                 Defaults to 'ALL'.
             forceNewFile(optional): Whether to create new output file.
             saveFormat(optional): Format for saved plot(s). Defaults to "eps".
             saveDPI(optional): DPI for saved plot(s). Defaults to 100.
             show(optional): Whether to show figure. Defaults to True.
             colorMap(optional): Pyplot colormap. Defaults to 'jet'.
             pyplotArgs(optional): Optional arguments for pyplot.
             figspecs(optional key args): Optional arguments for matplotlib
                 figure creation.
 
         Returns:
             Output filename and figure handle.
         """
         if isinstance(what, (str,)):
             if what.upper() == 'ALL':
                 what = ['ABS', 'PHASE', 'REAL', 'IMAG']
             else:
                 what = [what]
         what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'],
                          listname = self.name() + ".what", baselevel = 1)
         if len(what) == 0: return
         if 'figsize' not in figspecs.keys():
             figspecs['figsize'] = plt.figaspect(1. / len(what))
 
         if is_state or self.isCEye:
             idxs = np.arange(self.spacedim)
         else:
             idxs = np.arange(self.C.shape[0])
+        uP = u.reshape(idxs[-1] + 1, -1)
         if warping is not None:
             idxs = warping[0](np.arange(self.spacedim))
         subplotidx = 0
         fig = plt.figure(**figspecs)
         plt.set_cmap(colorMap)
         if 'ABS' in what:
             subplotidx = subplotidx + 1
             ax = fig.add_subplot(1, len(what), subplotidx)
-            ax.plot(idxs, np.abs(u).flatten(), **pyplotArgs)
+            ax.plot(idxs, np.abs(uP), **pyplotArgs)
             ax.set_title("|{0}|".format(name))
         if 'PHASE' in what:
             subplotidx = subplotidx + 1
             ax = fig.add_subplot(1, len(what), subplotidx)
-            ax.plot(idxs, np.angle(u).flatten(), **pyplotArgs)
+            ax.plot(idxs, np.angle(uP), **pyplotArgs)
             ax.set_title("phase({0})".format(name))
         if 'REAL' in what:
             subplotidx = subplotidx + 1
             ax = fig.add_subplot(1, len(what), subplotidx)
-            ax.plot(idxs, np.real(u).flatten(), **pyplotArgs)
+            ax.plot(idxs, np.real(uP), **pyplotArgs)
             ax.set_title("Re({0})".format(name))
         if 'IMAG' in what:
             subplotidx = subplotidx + 1
             ax = fig.add_subplot(1, len(what), subplotidx)
-            ax.plot(idxs, np.imag(u).flatten(), **pyplotArgs)
+            ax.plot(idxs, np.imag(uP), **pyplotArgs)
             ax.set_title("Im({0})".format(name))
         plt.tight_layout()
         if save is not None:
             save = save.strip()
             if forceNewFile:
                 fileOut = getNewFilename("{}_fig_".format(save), saveFormat)
             else:
                 fileOut = "{}_fig.{}".format(save, saveFormat)
             fig.savefig(fileOut, format = saveFormat, dpi = saveDPI)
         else: fileOut = None
         if show: plt.show()
         if fileOut is None: return fig
         return fig, fileOut
diff --git a/rrompy/hfengines/base/vector_fenics_engine_base.py b/rrompy/hfengines/base/vector_fenics_engine_base.py
index e0f7607..3b44ba9 100644
--- a/rrompy/hfengines/base/vector_fenics_engine_base.py
+++ b/rrompy/hfengines/base/vector_fenics_engine_base.py
@@ -1,148 +1,148 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import fenics as fen
 import numpy as np
 from matplotlib import pyplot as plt
 from .fenics_engine_base import FenicsEngineBase
 from rrompy.utilities.base.types import Np1D, List, strLst, Tuple, FigHandle
 from rrompy.utilities.base.data_structures import purgeList, getNewFilename
 from rrompy.solver.fenics import fenplot
 
 __all__ = ['VectorFenicsEngineBase']
 
 class VectorFenicsEngineBase(FenicsEngineBase):
     """Generic solver for parametric vector fenics problems."""
 
     def __init__(self, degree_threshold : int = np.inf,
                  verbosity : int = 10, timestamp : bool = True):
         super().__init__(degree_threshold = degree_threshold,
                          verbosity = verbosity, timestamp = timestamp)
         self.V = fen.VectorFunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
 
     def plot(self, u:Np1D, warping : List[callable] = None,
              is_state : bool = False, name : str = "u", save : str = None,
              what : strLst = 'all', forceNewFile : bool = True,
              saveFormat : str = "eps", saveDPI : int = 100,
              show : bool = True, colorMap : str = "jet",
              fenplotArgs : dict = {},
              **figspecs) -> Tuple[List[FigHandle], List[str]]:
         """
         Do some nice plots of the complex-valued function with given dofs.
 
         Args:
             u: numpy complex array with function dofs.
             warping(optional): Domain warping functions.
             is_state(optional): whether given u is value before multiplication
                 by c. Defaults to False.
             name(optional): Name to be shown as title of the plots. Defaults to
                 'u'.
             save(optional): Where to save plot(s). Defaults to None, i.e. no
                 saving.
             what(optional): Which plots to do. If list, can contain 'ABS',
                 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
                 Defaults to 'ALL'.
             forceNewFile(optional): Whether to create new output file.
             saveFormat(optional): Format for saved plot(s). Defaults to "eps".
             saveDPI(optional): DPI for saved plot(s). Defaults to.
             show(optional): Whether to show figure. Defaults to True.
             colorMap(optional): Pyplot colormap. Defaults to 'jet'.
             fenplotArgs(optional): Optional arguments for fenplot.
             figspecs(optional key args): Optional arguments for matplotlib
                 figure creation.
 
         Returns:
             List of output filenames and list of figure handles.
         """
         if not is_state and not self.isCEye:
             return super().plot(u, warping, False, name, save, what,
                                 forceNewFile, saveFormat, saveDPI, show,
-                                fenplotArgs, **figspecs)
+                                colorMap, fenplotArgs, **figspecs)
         if isinstance(what, (str,)):
             if what.upper() == 'ALL':
                 what = ['ABS', 'PHASE', 'REAL', 'IMAG']
             else:
                 what = [what]
         what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'],
                          listname = self.name() + ".what", baselevel = 1)
         if len(what) == 0: return
         if 'figsize' not in figspecs.keys():
             figspecs['figsize'] = plt.figaspect(1. / len(what))
 
         figs = [None] * self.V.num_sub_spaces()
         fileOut = [None] * self.V.num_sub_spaces()
         for j in range(self.V.num_sub_spaces()):
             II = self.V.sub(j).dofmap().dofs()
             Vj = self.V.sub(j).collapse()
             subplotidx = 0
             figs[j] = plt.figure(**figspecs)
             plt.set_cmap(colorMap)
             if 'ABS' in what:
                 uAb = fen.Function(Vj)
                 uAb.vector().set_local(np.abs(u[II]))
                 subplotidx = subplotidx + 1
                 ax = figs[j].add_subplot(1, len(what), subplotidx)
                 p = fenplot(uAb, warping = warping,
                             title = "|{}_comp{}|".format(name, j),
                             **fenplotArgs)
                 if self.V.mesh().geometric_dimension() > 1:
                     figs[j].colorbar(p, ax = ax)
             if 'PHASE' in what:
                 uPh = fen.Function(Vj)
                 uPh.vector().set_local(np.angle(u[II]))
                 subplotidx = subplotidx + 1
                 ax = figs[j].add_subplot(1, len(what), subplotidx)
                 p = fenplot(uPh, warping = warping,
                             title = "phase({}_comp{})".format(name, j),
                             **fenplotArgs)
                 if self.V.mesh().geometric_dimension() > 1:
                     figs[j].colorbar(p, ax = ax)
             if 'REAL' in what:
                 uRe = fen.Function(Vj)
                 uRe.vector().set_local(np.real(u[II]))
                 subplotidx = subplotidx + 1
                 ax = figs[j].add_subplot(1, len(what), subplotidx)
                 p = fenplot(uRe, warping = warping,
                             title = "Re({}_comp{})".format(name, j),
                             **fenplotArgs)
                 if self.V.mesh().geometric_dimension() > 1:
                     figs[j].colorbar(p, ax = ax)
             if 'IMAG' in what:
                 uIm = fen.Function(Vj)
                 uIm.vector().set_local(np.imag(u[II]))
                 subplotidx = subplotidx + 1
                 ax = figs[j].add_subplot(1, len(what), subplotidx)
                 p = fenplot(uIm, warping = warping,
                             title = "Im({}_comp{})".format(name, j),
                             **fenplotArgs)
                 if self.V.mesh().geometric_dimension() > 1:
                     figs[j].colorbar(p, ax = ax)
             plt.tight_layout()
             if save is not None:
                 save = save.strip()
                 if forceNewFile:
                     fileOut = getNewFilename(
                                           "{}_comp{}_fig_".format(save, j),
                                           saveFormat)
                 else:
                     fileOut = "{}_comp{}_fig.{}".format(save,j,saveFormat)
                 figs[j].savefig(fileOut, format = saveFormat,
                                       dpi = saveDPI)
             if show: plt.show()
         if fileOut[0] is None: return figs
         return figs, fileOut