Repository: sympy/sympy

Issue:
Plots are displaying incorrectly on the development documentation website
See: https://docs.sympy.org/dev/modules/physics/continuum_mechanics/beam.html#sympy.physics.continuum_mechanics.beam.Beam.plot_bending_moment

Other possible related issues:

- #18901
-  #19023
- #18974

Here are a couple of examples:

![image](https://user-images.githubusercontent.com/276007/78207512-c0455480-7456-11ea-94d2-6e2325a9fc70.png)

![image](https://user-images.githubusercontent.com/276007/78207555-de12b980-7456-11ea-9758-f726af7fc2c0.png)


Discussion:

Patch:
diff --git a/.travis.yml b/.travis.yml
index 346eaaadfd..f41dd8f3be 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -53,7 +53,7 @@ matrix:
       # all of the dependencies are supported on 3.8.
       env:
         - TEST_ASCII="true"
-        - TEST_OPT_DEPENDENCY="matchpy numpy scipy gmpy2 matplotlib<3.2 theano llvmlite autowrap cython wurlitzer python-symengine=0.5.1 tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle pyglet pycosat lfortran python-clang lxml"
+        - TEST_OPT_DEPENDENCY="matchpy numpy scipy gmpy2 matplotlib theano llvmlite autowrap cython wurlitzer python-symengine=0.5.1 tensorflow numexpr ipython antlr-python-runtime>=4.7,<4.8 antlr>=4.7,<4.8 cloudpickle pyglet pycosat lfortran python-clang lxml"
         - TEST_SAGE="true"
         - SYMPY_STRICT_COMPILER_CHECKS=1
       addons:
diff --git a/sympy/plotting/plot.py b/sympy/plotting/plot.py
index 9d9e95ed18..ae612e1442 100644
--- a/sympy/plotting/plot.py
+++ b/sympy/plotting/plot.py
@@ -144,29 +144,32 @@ class Plot(object):
     - surface_color : function which returns a float.
     """
 
-    def __init__(self, *args, **kwargs):
+    def __init__(self, *args,
+        title=None, xlabel=None, ylabel=None, aspect_ratio='auto',
+        xlim=None, ylim=None, axis_center='auto', axis=True,
+        xscale='linear', yscale='linear', legend=False, autoscale=True,
+        margin=0, annotations=None, markers=None, rectangles=None,
+        fill=None, backend='default', **kwargs):
         super(Plot, self).__init__()
 
         # Options for the graph as a whole.
         # The possible values for each option are described in the docstring of
         # Plot. They are based purely on convention, no checking is done.
-        self.title = None
-        self.xlabel = None
-        self.ylabel = None
-        self.aspect_ratio = 'auto'
-        self.xlim = None
-        self.ylim = None
-        self.axis_center = 'auto'
-        self.axis = True
-        self.xscale = 'linear'
-        self.yscale = 'linear'
-        self.legend = False
-        self.autoscale = True
-        self.margin = 0
-        self.annotations = None
-        self.markers = None
-        self.rectangles = None
-        self.fill = None
+        self.title = title
+        self.xlabel = xlabel
+        self.ylabel = ylabel
+        self.aspect_ratio = aspect_ratio
+        self.axis_center = axis_center
+        self.axis = axis
+        self.xscale = xscale
+        self.yscale = yscale
+        self.legend = legend
+        self.autoscale = autoscale
+        self.margin = margin
+        self.annotations = annotations
+        self.markers = markers
+        self.rectangles = rectangles
+        self.fill = fill
 
         # Contains the data objects to be plotted. The backend should be smart
         # enough to iterate over this list.
@@ -176,13 +179,32 @@ def __init__(self, *args, **kwargs):
         # The backend type. On every show() a new backend instance is created
         # in self._backend which is tightly coupled to the Plot instance
         # (thanks to the parent attribute of the backend).
-        self.backend = plot_backends[kwargs.pop('backend', 'default')]
+        self.backend = plot_backends[backend]
 
+        is_real = \
+            lambda lim: all(getattr(i, 'is_real', True) for i in lim)
+        is_finite = \
+            lambda lim: all(getattr(i, 'is_finite', True) for i in lim)
+
+        self.xlim = None
+        self.ylim = None
+        if xlim:
+            if not is_real(xlim):
+                raise ValueError(
+                "All numbers from xlim={} must be real".format(xlim))
+            if not is_finite(xlim):
+                raise ValueError(
+                "All numbers from xlim={} must be finite".format(xlim))
+            self.xlim = (float(xlim[0]), float(xlim[1]))
+        if ylim:
+            if not is_real(ylim):
+                raise ValueError(
+                "All numbers from ylim={} must be real".format(ylim))
+            if not is_finite(ylim):
+                raise ValueError(
+                "All numbers from ylim={} must be finite".format(ylim))
+            self.ylim = (float(ylim[0]), float(ylim[1]))
 
-        # The keyword arguments should only contain options for the plot.
-        for key, val in kwargs.items():
-            if hasattr(self, key):
-                setattr(self, key, val)
 
     def show(self):
         # TODO move this to the backend (also for save)
@@ -847,14 +869,28 @@ def get_parameter_points(self):
         return np.linspace(self.start, self.end, num=self.nb_of_points)
 
     def get_points(self):
+        np = import_module('numpy')
         param = self.get_parameter_points()
         fx = vectorized_lambdify([self.var], self.expr_x)
         fy = vectorized_lambdify([self.var], self.expr_y)
         fz = vectorized_lambdify([self.var], self.expr_z)
+
         list_x = fx(param)
         list_y = fy(param)
         list_z = fz(param)
-        return (list_x, list_y, list_z)
+
+        list_x = np.array(list_x, dtype=np.float64)
+        list_y = np.array(list_y, dtype=np.float64)
+        list_z = np.array(list_z, dtype=np.float64)
+
+        list_x = np.ma.masked_invalid(list_x)
+        list_y = np.ma.masked_invalid(list_y)
+        list_z = np.ma.masked_invalid(list_z)
+
+        self._xlim = (np.amin(list_x), np.amax(list_x))
+        self._ylim = (np.amin(list_y), np.amax(list_y))
+        self._zlim = (np.amin(list_z), np.amax(list_z))
+        return list_x, list_y, list_z
 
 
 ### Surfaces
@@ -906,6 +942,9 @@ def __init__(self, expr, var_start_end_x, var_start_end_y, **kwargs):
         self.nb_of_points_y = kwargs.get('nb_of_points_y', 50)
         self.surface_color = kwargs.get('surface_color', None)
 
+        self._xlim = (self.start_x, self.end_x)
+        self._ylim = (self.start_y, self.end_y)
+
     def __str__(self):
         return ('cartesian surface: %s for'
                 ' %s over %s and %s over %s') % (
@@ -922,7 +961,11 @@ def get_meshes(self):
                                      np.linspace(self.start_y, self.end_y,
                                                  num=self.nb_of_points_y))
         f = vectorized_lambdify((self.var_x, self.var_y), self.expr)
-        return (mesh_x, mesh_y, f(mesh_x, mesh_y))
+        mesh_z = f(mesh_x, mesh_y)
+        mesh_z = np.array(mesh_z, dtype=np.float64)
+        mesh_z = np.ma.masked_invalid(mesh_z)
+        self._zlim = (np.amin(mesh_z), np.amax(mesh_z))
+        return mesh_x, mesh_y, mesh_z
 
 
 class ParametricSurfaceSeries(SurfaceBaseSeries):
@@ -967,11 +1010,30 @@ def get_parameter_meshes(self):
                                        num=self.nb_of_points_v))
 
     def get_meshes(self):
+        np = import_module('numpy')
+
         mesh_u, mesh_v = self.get_parameter_meshes()
         fx = vectorized_lambdify((self.var_u, self.var_v), self.expr_x)
         fy = vectorized_lambdify((self.var_u, self.var_v), self.expr_y)
         fz = vectorized_lambdify((self.var_u, self.var_v), self.expr_z)
-        return (fx(mesh_u, mesh_v), fy(mesh_u, mesh_v), fz(mesh_u, mesh_v))
+
+        mesh_x = fx(mesh_u, mesh_v)
+        mesh_y = fy(mesh_u, mesh_v)
+        mesh_z = fz(mesh_u, mesh_v)
+
+        mesh_x = np.array(mesh_x, dtype=np.float64)
+        mesh_y = np.array(mesh_y, dtype=np.float64)
+        mesh_z = np.array(mesh_z, dtype=np.float64)
+
+        mesh_x = np.ma.masked_invalid(mesh_x)
+        mesh_y = np.ma.masked_invalid(mesh_y)
+        mesh_z = np.ma.masked_invalid(mesh_z)
+
+        self._xlim = (np.amin(mesh_x), np.amax(mesh_x))
+        self._ylim = (np.amin(mesh_y), np.amax(mesh_y))
+        self._zlim = (np.amin(mesh_z), np.amax(mesh_z))
+
+        return mesh_x, mesh_y, mesh_z
 
 
 ### Contours
@@ -996,6 +1058,9 @@ def __init__(self, expr, var_start_end_x, var_start_end_y):
 
         self.get_points = self.get_meshes
 
+        self._xlim = (self.start_x, self.end_x)
+        self._ylim = (self.start_y, self.end_y)
+
     def __str__(self):
         return ('contour: %s for '
                 '%s over %s and %s over %s') % (
@@ -1068,12 +1133,18 @@ def __init__(self, parent):
                 self.ax[i].spines['right'].set_color('none')
                 self.ax[i].spines['bottom'].set_position('zero')
                 self.ax[i].spines['top'].set_color('none')
-                self.ax[i].spines['left'].set_smart_bounds(True)
-                self.ax[i].spines['bottom'].set_smart_bounds(False)
                 self.ax[i].xaxis.set_ticks_position('bottom')
                 self.ax[i].yaxis.set_ticks_position('left')
 
     def _process_series(self, series, ax, parent):
+        np = import_module('numpy')
+        mpl_toolkits = import_module(
+            'mpl_toolkits', import_kwargs={'fromlist': ['mplot3d']})
+
+        # XXX Workaround for matplotlib issue
+        # https://github.com/matplotlib/matplotlib/issues/17130
+        xlims, ylims, zlims = [], [], []
+
         for s in series:
             # Create the collections
             if s.is_2Dline:
@@ -1083,24 +1154,22 @@ def _process_series(self, series, ax, parent):
                 ax.contour(*s.get_meshes())
             elif s.is_3Dline:
                 # TODO too complicated, I blame matplotlib
-                mpl_toolkits = import_module('mpl_toolkits',
-                    import_kwargs={'fromlist': ['mplot3d']})
                 art3d = mpl_toolkits.mplot3d.art3d
                 collection = art3d.Line3DCollection(s.get_segments())
                 ax.add_collection(collection)
                 x, y, z = s.get_points()
-                ax.set_xlim((min(x), max(x)))
-                ax.set_ylim((min(y), max(y)))
-                ax.set_zlim((min(z), max(z)))
+                xlims.append(s._xlim)
+                ylims.append(s._ylim)
+                zlims.append(s._zlim)
             elif s.is_3Dsurface:
                 x, y, z = s.get_meshes()
                 collection = ax.plot_surface(x, y, z,
                     cmap=getattr(self.cm, 'viridis', self.cm.jet),
                     rstride=1, cstride=1, linewidth=0.1)
+                xlims.append(s._xlim)
+                ylims.append(s._ylim)
+                zlims.append(s._zlim)
             elif s.is_implicit:
-                # Smart bounds have to be set to False for implicit plots.
-                ax.spines['left'].set_smart_bounds(False)
-                ax.spines['bottom'].set_smart_bounds(False)
                 points = s.get_raster()
                 if len(points) == 2:
                     # interval math plotting
@@ -1118,9 +1187,10 @@ def _process_series(self, series, ax, parent):
                     else:
                         ax.contourf(xarray, yarray, zarray, cmap=colormap)
             else:
-                raise ValueError('The matplotlib backend supports only '
-                                 'is_2Dline, is_3Dline, is_3Dsurface and '
-                                 'is_contour objects.')
+                raise NotImplementedError(
+                    '{} is not supported in the sympy plotting module '
+                    'with matplotlib backend. Please report this issue.'
+                    .format(ax))
 
             # Customise the collections with the corresponding per-series
             # options.
@@ -1142,12 +1212,38 @@ def _process_series(self, series, ax, parent):
                 else:
                     collection.set_color(s.surface_color)
 
+        Axes3D = mpl_toolkits.mplot3d.Axes3D
+        if not isinstance(ax, Axes3D):
+            ax.autoscale_view(
+                scalex=ax.get_autoscalex_on(),
+                scaley=ax.get_autoscaley_on())
+        else:
+            # XXX Workaround for matplotlib issue
+            # https://github.com/matplotlib/matplotlib/issues/17130
+            if xlims:
+                xlims = np.array(xlims)
+                xlim = (np.amin(xlims[:, 0]), np.amax(xlims[:, 1]))
+                ax.set_xlim(xlim)
+            else:
+                ax.set_xlim([0, 1])
+
+            if ylims:
+                ylims = np.array(ylims)
+                ylim = (np.amin(ylims[:, 0]), np.amax(ylims[:, 1]))
+                ax.set_ylim(ylim)
+            else:
+                ax.set_ylim([0, 1])
+
+            if zlims:
+                zlims = np.array(zlims)
+                zlim = (np.amin(zlims[:, 0]), np.amax(zlims[:, 1]))
+                ax.set_zlim(zlim)
+            else:
+                ax.set_zlim([0, 1])
+
         # Set global options.
         # TODO The 3D stuff
         # XXX The order of those is important.
-        mpl_toolkits = import_module('mpl_toolkits',
-            import_kwargs={'fromlist': ['mplot3d']})
-        Axes3D = mpl_toolkits.mplot3d.Axes3D
         if parent.xscale and not isinstance(ax, Axes3D):
             ax.set_xscale(parent.xscale)
         if parent.yscale and not isinstance(ax, Axes3D):
@@ -1205,38 +1301,9 @@ def _process_series(self, series, ax, parent):
         # xlim and ylim shoulld always be set at last so that plot limits
         # doesn't get altered during the process.
         if parent.xlim:
-            from sympy.core.basic import Basic
-            xlim = parent.xlim
-            if any(isinstance(i, Basic) and not i.is_real for i in xlim):
-                raise ValueError(
-                "All numbers from xlim={} must be real".format(xlim))
-            if any(isinstance(i, Basic) and not i.is_finite for i in xlim):
-                raise ValueError(
-                "All numbers from xlim={} must be finite".format(xlim))
-            xlim = (float(i) for i in xlim)
-            ax.set_xlim(xlim)
-        else:
-            if parent._series and all(isinstance(s, LineOver1DRangeSeries) for s in parent._series):
-                starts = [s.start for s in parent._series]
-                ends = [s.end for s in parent._series]
-                ax.set_xlim(min(starts), max(ends))
-
+            ax.set_xlim(parent.xlim)
         if parent.ylim:
-            from sympy.core.basic import Basic
-            ylim = parent.ylim
-            if any(isinstance(i,Basic) and not i.is_real for i in ylim):
-                raise ValueError(
-                "All numbers from ylim={} must be real".format(ylim))
-            if any(isinstance(i,Basic) and not i.is_finite for i in ylim):
-                raise ValueError(
-                "All numbers from ylim={} must be finite".format(ylim))
-            ylim = (float(i) for i in ylim)
-            ax.set_ylim(ylim)
-
-        if not isinstance(ax, Axes3D):
-            ax.autoscale_view(
-                scalex=ax.get_autoscalex_on(),
-                scaley=ax.get_autoscaley_on())
+            ax.set_ylim(parent.ylim)
 
 
     def process_series(self):


Test patch:
diff --git a/sympy/plotting/tests/test_experimental_lambdify.py b/sympy/plotting/tests/test_experimental_lambdify.py
index ea59a437d5..95839d6687 100644
--- a/sympy/plotting/tests/test_experimental_lambdify.py
+++ b/sympy/plotting/tests/test_experimental_lambdify.py
@@ -1,9 +1,25 @@
-from sympy.core.symbol import symbols
+from sympy.core.symbol import symbols, Symbol
+from sympy.functions import Max
 from sympy.plotting.experimental_lambdify import experimental_lambdify
 from sympy.plotting.intervalmath.interval_arithmetic import \
     interval, intervalMembership
 
 
+# Tests for exception handling in experimental_lambdify
+def test_experimental_lambify():
+    x = Symbol('x')
+    f = experimental_lambdify([x], Max(x, 5))
+    # XXX should f be tested? If f(2) is attempted, an
+    # error is raised because a complex produced during wrapping of the arg
+    # is being compared with an int.
+    assert Max(2, 5) == 5
+    assert Max(5, 7) == 7
+
+    x = Symbol('x-3')
+    f = experimental_lambdify([x], x + 1)
+    assert f(1) == 2
+
+
 def test_composite_boolean_region():
     x, y = symbols('x y')
 
diff --git a/sympy/plotting/tests/test_plot.py b/sympy/plotting/tests/test_plot.py
index 48ea9db2a4..ef6bf678f3 100644
--- a/sympy/plotting/tests/test_plot.py
+++ b/sympy/plotting/tests/test_plot.py
@@ -1,271 +1,314 @@
-from typing import List
-
-from sympy import (pi, sin, cos, Symbol, Integral, Sum, sqrt, log, exp, Ne,
-                   oo, LambertW, I, meijerg, exp_polar, Max, Piecewise, And,
-                   real_root)
-from sympy.plotting import (plot, plot_parametric, plot3d_parametric_line,
-                            plot3d, plot3d_parametric_surface)
-from sympy.plotting.plot import (unset_show, plot_contour, PlotGrid,
-                            DefaultBackend, MatplotlibBackend, TextBackend)
-from sympy.utilities import lambdify as lambdify_
-from sympy.testing.pytest import skip, raises, warns
-from sympy.plotting.experimental_lambdify import lambdify
-from sympy.external import import_module
-
-from tempfile import NamedTemporaryFile
 import os
+from tempfile import TemporaryDirectory
 
-unset_show()
+from sympy import (
+    pi, sin, cos, Symbol, Integral, Sum, sqrt, log, exp, Ne, oo, LambertW, I,
+    meijerg, exp_polar, Piecewise, And, real_root)
+from sympy.core.singleton import S
+from sympy.core.sympify import sympify
+from sympy.external import import_module
+from sympy.plotting.plot import (
+    Plot, plot, plot_parametric, plot3d_parametric_line, plot3d,
+    plot3d_parametric_surface)
+from sympy.plotting.plot import (
+    unset_show, plot_contour, PlotGrid, DefaultBackend, MatplotlibBackend,
+    TextBackend)
+from sympy.testing.pytest import skip, raises, warns
+from sympy.utilities import lambdify as lambdify_
 
 
-# XXX: We could implement this as a context manager instead
-# That would need rewriting the plot_and_save() function
-# entirely
-class TmpFileManager:
+unset_show()
 
-    tmp_files = []  # type: List[str]
 
-    @classmethod
-    def tmp_file(cls, name=''):
-        cls.tmp_files.append(NamedTemporaryFile(prefix=name, suffix='.png').name)
-        return cls.tmp_files[-1]
+matplotlib = import_module(
+    'matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
 
-    @classmethod
-    def cleanup(cls):
-        for file in cls.tmp_files:
-            try:
-                os.remove(file)
-            except OSError:
-                # If the file doesn't exist, for instance, if the test failed.
-                pass
 
-def plot_and_save_1(name):
-    tmp_file = TmpFileManager.tmp_file
+def test_plot_and_save_1():
+    if not matplotlib:
+        skip("Matplotlib not the default backend")
 
     x = Symbol('x')
     y = Symbol('y')
 
-    ###
-    # Examples from the 'introduction' notebook
-    ###
+    with TemporaryDirectory(prefix='sympy_') as tmpdir:
+        ###
+        # Examples from the 'introduction' notebook
+        ###
+        p = plot(x, legend=True, label='f1')
+        p = plot(x*sin(x), x*cos(x), label='f2')
+        p.extend(p)
+        p[0].line_color = lambda a: a
+        p[1].line_color = 'b'
+        p.title = 'Big title'
+        p.xlabel = 'the x axis'
+        p[1].label = 'straight line'
+        p.legend = True
+        p.aspect_ratio = (1, 1)
+        p.xlim = (-15, 20)
+        filename = 'test_basic_options_and_colors.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot(x, legend=True, label='f1')
-    p = plot(x*sin(x), x*cos(x), label='f2')
-    p.extend(p)
-    p[0].line_color = lambda a: a
-    p[1].line_color = 'b'
-    p.title = 'Big title'
-    p.xlabel = 'the x axis'
-    p[1].label = 'straight line'
-    p.legend = True
-    p.aspect_ratio = (1, 1)
-    p.xlim = (-15, 20)
-    p.save(tmp_file('%s_basic_options_and_colors' % name))
-    p._backend.close()
+        p.extend(plot(x + 1))
+        p.append(plot(x + 3, x**2)[1])
+        filename = 'test_plot_extend_append.png'
+        p.save(os.path.join(tmpdir, filename))
 
-    p.extend(plot(x + 1))
-    p.append(plot(x + 3, x**2)[1])
-    p.save(tmp_file('%s_plot_extend_append' % name))
+        p[2] = plot(x**2, (x, -2, 3))
+        filename = 'test_plot_setitem.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p[2] = plot(x**2, (x, -2, 3))
-    p.save(tmp_file('%s_plot_setitem' % name))
-    p._backend.close()
+        p = plot(sin(x), (x, -2*pi, 4*pi))
+        filename = 'test_line_explicit.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot(sin(x), (x, -2*pi, 4*pi))
-    p.save(tmp_file('%s_line_explicit' % name))
-    p._backend.close()
+        p = plot(sin(x))
+        filename = 'test_line_default_range.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot(sin(x))
-    p.save(tmp_file('%s_line_default_range' % name))
-    p._backend.close()
+        p = plot((x**2, (x, -5, 5)), (x**3, (x, -3, 3)))
+        filename = 'test_line_multiple_range.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot((x**2, (x, -5, 5)), (x**3, (x, -3, 3)))
-    p.save(tmp_file('%s_line_multiple_range' % name))
-    p._backend.close()
+        raises(ValueError, lambda: plot(x, y))
 
-    raises(ValueError, lambda: plot(x, y))
+        #Piecewise plots
+        p = plot(Piecewise((1, x > 0), (0, True)), (x, -1, 1))
+        filename = 'test_plot_piecewise.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    #Piecewise plots
-    p = plot(Piecewise((1, x > 0), (0, True)), (x, -1, 1))
-    p.save(tmp_file('%s_plot_piecewise' % name))
-    p._backend.close()
+        p = plot(Piecewise((x, x < 1), (x**2, True)), (x, -3, 3))
+        filename = 'test_plot_piecewise_2.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot(Piecewise((x, x < 1), (x**2, True)), (x, -3, 3))
-    p.save(tmp_file('%s_plot_piecewise_2' % name))
-    p._backend.close()
+        # test issue 7471
+        p1 = plot(x)
+        p2 = plot(3)
+        p1.extend(p2)
+        filename = 'test_horizontal_line.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # test issue 7471
-    p1 = plot(x)
-    p2 = plot(3)
-    p1.extend(p2)
-    p.save(tmp_file('%s_horizontal_line' % name))
-    p._backend.close()
+        # test issue 10925
+        f = Piecewise((-1, x < -1), (x, And(-1 <= x, x < 0)), \
+            (x**2, And(0 <= x, x < 1)), (x**3, x >= 1))
+        p = plot(f, (x, -3, 3))
+        filename = 'test_plot_piecewise_3.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # test issue 10925
-    f = Piecewise((-1, x < -1), (x, And(-1 <= x, x < 0)), \
-        (x**2, And(0 <= x, x < 1)), (x**3, x >= 1))
-    p = plot(f, (x, -3, 3))
-    p.save(tmp_file('%s_plot_piecewise_3' % name))
-    p._backend.close()
 
-def plot_and_save_2(name):
-    tmp_file = TmpFileManager.tmp_file
+def test_plot_and_save_2():
+    if not matplotlib:
+        skip("Matplotlib not the default backend")
 
     x = Symbol('x')
     y = Symbol('y')
     z = Symbol('z')
 
-    #parametric 2d plots.
-    #Single plot with default range.
-    plot_parametric(sin(x), cos(x)).save(tmp_file())
+    with TemporaryDirectory(prefix='sympy_') as tmpdir:
+        #parametric 2d plots.
+        #Single plot with default range.
+        p = plot_parametric(sin(x), cos(x))
+        filename = 'test_parametric.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    #Single plot with range.
-    p = plot_parametric(sin(x), cos(x), (x, -5, 5), legend=True, label='parametric_plot')
-    p.save(tmp_file('%s_parametric_range' % name))
-    p._backend.close()
+        #Single plot with range.
+        p = plot_parametric(
+            sin(x), cos(x), (x, -5, 5), legend=True, label='parametric_plot')
+        filename = 'test_parametric_range.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    #Multiple plots with same range.
-    p = plot_parametric((sin(x), cos(x)), (x, sin(x)))
-    p.save(tmp_file('%s_parametric_multiple' % name))
-    p._backend.close()
+        #Multiple plots with same range.
+        p = plot_parametric((sin(x), cos(x)), (x, sin(x)))
+        filename = 'test_parametric_multiple.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    #Multiple plots with different ranges.
-    p = plot_parametric((sin(x), cos(x), (x, -3, 3)), (x, sin(x), (x, -5, 5)))
-    p.save(tmp_file('%s_parametric_multiple_ranges' % name))
-    p._backend.close()
+        #Multiple plots with different ranges.
+        p = plot_parametric(
+            (sin(x), cos(x), (x, -3, 3)), (x, sin(x), (x, -5, 5)))
+        filename = 'test_parametric_multiple_ranges.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    #depth of recursion specified.
-    p = plot_parametric(x, sin(x), depth=13)
-    p.save(tmp_file('%s_recursion_depth' % name))
-    p._backend.close()
+        #depth of recursion specified.
+        p = plot_parametric(x, sin(x), depth=13)
+        filename = 'test_recursion_depth.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    #No adaptive sampling.
-    p = plot_parametric(cos(x), sin(x), adaptive=False, nb_of_points=500)
-    p.save(tmp_file('%s_adaptive' % name))
-    p._backend.close()
+        #No adaptive sampling.
+        p = plot_parametric(cos(x), sin(x), adaptive=False, nb_of_points=500)
+        filename = 'test_adaptive.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    #3d parametric plots
-    p = plot3d_parametric_line(sin(x), cos(x), x, legend=True, label='3d_parametric_plot')
-    p.save(tmp_file('%s_3d_line' % name))
-    p._backend.close()
+        #3d parametric plots
+        p = plot3d_parametric_line(
+            sin(x), cos(x), x, legend=True, label='3d_parametric_plot')
+        filename = 'test_3d_line.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot3d_parametric_line(
-        (sin(x), cos(x), x, (x, -5, 5)), (cos(x), sin(x), x, (x, -3, 3)))
-    p.save(tmp_file('%s_3d_line_multiple' % name))
-    p._backend.close()
+        p = plot3d_parametric_line(
+            (sin(x), cos(x), x, (x, -5, 5)), (cos(x), sin(x), x, (x, -3, 3)))
+        filename = 'test_3d_line_multiple.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot3d_parametric_line(sin(x), cos(x), x, nb_of_points=30)
-    p.save(tmp_file('%s_3d_line_points' % name))
-    p._backend.close()
+        p = plot3d_parametric_line(sin(x), cos(x), x, nb_of_points=30)
+        filename = 'test_3d_line_points.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # 3d surface single plot.
-    p = plot3d(x * y)
-    p.save(tmp_file('%s_surface' % name))
-    p._backend.close()
+        # 3d surface single plot.
+        p = plot3d(x * y)
+        filename = 'test_surface.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # Multiple 3D plots with same range.
-    p = plot3d(-x * y, x * y, (x, -5, 5))
-    p.save(tmp_file('%s_surface_multiple' % name))
-    p._backend.close()
+        # Multiple 3D plots with same range.
+        p = plot3d(-x * y, x * y, (x, -5, 5))
+        filename = 'test_surface_multiple.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # Multiple 3D plots with different ranges.
-    p = plot3d(
-        (x * y, (x, -3, 3), (y, -3, 3)), (-x * y, (x, -3, 3), (y, -3, 3)))
-    p.save(tmp_file('%s_surface_multiple_ranges' % name))
-    p._backend.close()
+        # Multiple 3D plots with different ranges.
+        p = plot3d(
+            (x * y, (x, -3, 3), (y, -3, 3)), (-x * y, (x, -3, 3), (y, -3, 3)))
+        filename = 'test_surface_multiple_ranges.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # Single Parametric 3D plot
-    p = plot3d_parametric_surface(sin(x + y), cos(x - y), x - y)
-    p.save(tmp_file('%s_parametric_surface' % name))
-    p._backend.close()
+        # Single Parametric 3D plot
+        p = plot3d_parametric_surface(sin(x + y), cos(x - y), x - y)
+        filename = 'test_parametric_surface.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # Multiple Parametric 3D plots.
-    p = plot3d_parametric_surface(
-        (x*sin(z), x*cos(z), z, (x, -5, 5), (z, -5, 5)),
-        (sin(x + y), cos(x - y), x - y, (x, -5, 5), (y, -5, 5)))
-    p.save(tmp_file('%s_parametric_surface' % name))
-    p._backend.close()
+        # Multiple Parametric 3D plots.
+        p = plot3d_parametric_surface(
+            (x*sin(z), x*cos(z), z, (x, -5, 5), (z, -5, 5)),
+            (sin(x + y), cos(x - y), x - y, (x, -5, 5), (y, -5, 5)))
+        filename = 'test_parametric_surface.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # Single Contour plot.
-    p = plot_contour(sin(x)*sin(y), (x, -5, 5), (y, -5, 5))
-    p.save(tmp_file('%s_contour_plot' % name))
-    p._backend.close()
+        # Single Contour plot.
+        p = plot_contour(sin(x)*sin(y), (x, -5, 5), (y, -5, 5))
+        filename = 'test_contour_plot.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # Multiple Contour plots with same range.
-    p = plot_contour(x**2 + y**2, x**3 + y**3, (x, -5, 5), (y, -5, 5))
-    p.save(tmp_file('%s_contour_plot' % name))
-    p._backend.close()
+        # Multiple Contour plots with same range.
+        p = plot_contour(x**2 + y**2, x**3 + y**3, (x, -5, 5), (y, -5, 5))
+        filename = 'test_contour_plot.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
+
+        # Multiple Contour plots with different range.
+        p = plot_contour(
+            (x**2 + y**2, (x, -5, 5), (y, -5, 5)),
+            (x**3 + y**3, (x, -3, 3), (y, -3, 3)))
+        filename = 'test_contour_plot.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    # Multiple Contour plots with different range.
-    p = plot_contour((x**2 + y**2, (x, -5, 5), (y, -5, 5)), (x**3 + y**3, (x, -3, 3), (y, -3, 3)))
-    p.save(tmp_file('%s_contour_plot' % name))
-    p._backend.close()
 
-def plot_and_save_3(name):
-    tmp_file = TmpFileManager.tmp_file
+def test_plot_and_save_3():
+    if not matplotlib:
+        skip("Matplotlib not the default backend")
 
     x = Symbol('x')
     y = Symbol('y')
     z = Symbol('z')
 
-    ###
-    # Examples from the 'colors' notebook
-    ###
+    with TemporaryDirectory(prefix='sympy_') as tmpdir:
+        ###
+        # Examples from the 'colors' notebook
+        ###
 
-    p = plot(sin(x))
-    p[0].line_color = lambda a: a
-    p.save(tmp_file('%s_colors_line_arity1' % name))
+        p = plot(sin(x))
+        p[0].line_color = lambda a: a
+        filename = 'test_colors_line_arity1.png'
+        p.save(os.path.join(tmpdir, filename))
 
-    p[0].line_color = lambda a, b: b
-    p.save(tmp_file('%s_colors_line_arity2' % name))
-    p._backend.close()
+        p[0].line_color = lambda a, b: b
+        filename = 'test_colors_line_arity2.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot(x*sin(x), x*cos(x), (x, 0, 10))
-    p[0].line_color = lambda a: a
-    p.save(tmp_file('%s_colors_param_line_arity1' % name))
+        p = plot(x*sin(x), x*cos(x), (x, 0, 10))
+        p[0].line_color = lambda a: a
+        filename = 'test_colors_param_line_arity1.png'
+        p.save(os.path.join(tmpdir, filename))
 
-    p[0].line_color = lambda a, b: a
-    p.save(tmp_file('%s_colors_param_line_arity2a' % name))
+        p[0].line_color = lambda a, b: a
+        filename = 'test_colors_param_line_arity1.png'
+        p.save(os.path.join(tmpdir, filename))
 
-    p[0].line_color = lambda a, b: b
-    p.save(tmp_file('%s_colors_param_line_arity2b' % name))
-    p._backend.close()
+        p[0].line_color = lambda a, b: b
+        filename = 'test_colors_param_line_arity2b.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot3d_parametric_line(sin(x) + 0.1*sin(x)*cos(7*x),
-             cos(x) + 0.1*cos(x)*cos(7*x),
-        0.1*sin(7*x),
-        (x, 0, 2*pi))
-    p[0].line_color = lambdify_(x, sin(4*x))
-    p.save(tmp_file('%s_colors_3d_line_arity1' % name))
-    p[0].line_color = lambda a, b: b
-    p.save(tmp_file('%s_colors_3d_line_arity2' % name))
-    p[0].line_color = lambda a, b, c: c
-    p.save(tmp_file('%s_colors_3d_line_arity3' % name))
-    p._backend.close()
+        p = plot3d_parametric_line(sin(x) + 0.1*sin(x)*cos(7*x),
+                cos(x) + 0.1*cos(x)*cos(7*x),
+            0.1*sin(7*x),
+            (x, 0, 2*pi))
+        p[0].line_color = lambdify_(x, sin(4*x))
+        filename = 'test_colors_3d_line_arity1.png'
+        p.save(os.path.join(tmpdir, filename))
+        p[0].line_color = lambda a, b: b
+        filename = 'test_colors_3d_line_arity2.png'
+        p.save(os.path.join(tmpdir, filename))
+        p[0].line_color = lambda a, b, c: c
+        filename = 'test_colors_3d_line_arity3.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot3d(sin(x)*y, (x, 0, 6*pi), (y, -5, 5))
-    p[0].surface_color = lambda a: a
-    p.save(tmp_file('%s_colors_surface_arity1' % name))
-    p[0].surface_color = lambda a, b: b
-    p.save(tmp_file('%s_colors_surface_arity2' % name))
-    p[0].surface_color = lambda a, b, c: c
-    p.save(tmp_file('%s_colors_surface_arity3a' % name))
-    p[0].surface_color = lambdify_((x, y, z), sqrt((x - 3*pi)**2 + y**2))
-    p.save(tmp_file('%s_colors_surface_arity3b' % name))
-    p._backend.close()
+        p = plot3d(sin(x)*y, (x, 0, 6*pi), (y, -5, 5))
+        p[0].surface_color = lambda a: a
+        filename = 'test_colors_surface_arity1.png'
+        p.save(os.path.join(tmpdir, filename))
+        p[0].surface_color = lambda a, b: b
+        filename = 'test_colors_surface_arity2.png'
+        p.save(os.path.join(tmpdir, filename))
+        p[0].surface_color = lambda a, b, c: c
+        filename = 'test_colors_surface_arity3a.png'
+        p.save(os.path.join(tmpdir, filename))
+        p[0].surface_color = lambdify_((x, y, z), sqrt((x - 3*pi)**2 + y**2))
+        filename = 'test_colors_surface_arity3b.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    p = plot3d_parametric_surface(x * cos(4 * y), x * sin(4 * y), y,
-             (x, -1, 1), (y, -1, 1))
-    p[0].surface_color = lambda a: a
-    p.save(tmp_file('%s_colors_param_surf_arity1' % name))
-    p[0].surface_color = lambda a, b: a*b
-    p.save(tmp_file('%s_colors_param_surf_arity2' % name))
-    p[0].surface_color = lambdify_((x, y, z), sqrt(x**2 + y**2 + z**2))
-    p.save(tmp_file('%s_colors_param_surf_arity3' % name))
-    p._backend.close()
+        p = plot3d_parametric_surface(x * cos(4 * y), x * sin(4 * y), y,
+                (x, -1, 1), (y, -1, 1))
+        p[0].surface_color = lambda a: a
+        filename = 'test_colors_param_surf_arity1.png'
+        p.save(os.path.join(tmpdir, filename))
+        p[0].surface_color = lambda a, b: a*b
+        filename = 'test_colors_param_surf_arity2.png'
+        p.save(os.path.join(tmpdir, filename))
+        p[0].surface_color = lambdify_((x, y, z), sqrt(x**2 + y**2 + z**2))
+        filename = 'test_colors_param_surf_arity3.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-def plot_and_save_4(name):
-    tmp_file = TmpFileManager.tmp_file
+
+def test_plot_and_save_4():
+    if not matplotlib:
+        skip("Matplotlib not the default backend")
 
     x = Symbol('x')
     y = Symbol('y')
@@ -279,183 +322,109 @@ def plot_and_save_4(name):
     # report this as a bug." It has to use the fallback because using evalf()
     # is the only way to evaluate the integral. We should perhaps just remove
     # that warning.
-
-    with warns(UserWarning, match="The evaluation of the expression is problematic"):
-        i = Integral(log((sin(x)**2 + 1)*sqrt(x**2 + 1)), (x, 0, y))
-        p = plot(i, (y, 1, 5))
-        p.save(tmp_file('%s_advanced_integral' % name))
-        p._backend.close()
-
-def plot_and_save_5(name):
-    tmp_file = TmpFileManager.tmp_file
-
-    x = Symbol('x')
-    y = Symbol('y')
-
-    s = Sum(1/x**y, (x, 1, oo))
-    p = plot(s, (y, 2, 10))
-    p.save(tmp_file('%s_advanced_inf_sum' % name))
-    p._backend.close()
-
-    p = plot(Sum(1/x, (x, 1, y)), (y, 2, 10), show=False)
-    p[0].only_integers = True
-    p[0].steps = True
-    p.save(tmp_file('%s_advanced_fin_sum' % name))
-    p._backend.close()
-
-def plot_and_save_6(name):
-    tmp_file = TmpFileManager.tmp_file
-
-    x = Symbol('x')
-
-    ###
-    # Test expressions that can not be translated to np and generate complex
-    # results.
-    ###
-    plot(sin(x) + I*cos(x)).save(tmp_file())
-    plot(sqrt(sqrt(-x))).save(tmp_file())
-    plot(LambertW(x)).save(tmp_file())
-    plot(sqrt(LambertW(x))).save(tmp_file())
-
-    #Characteristic function of a StudentT distribution with nu=10
-    plot((meijerg(((1 / 2,), ()), ((5, 0, 1 / 2), ()), 5 * x**2 * exp_polar(-I*pi)/2)
-            + meijerg(((1/2,), ()), ((5, 0, 1/2), ()),
-                5*x**2 * exp_polar(I*pi)/2)) / (48 * pi), (x, 1e-6, 1e-2)).save(tmp_file())
-
-
-def plotgrid_and_save(name):
-    tmp_file = TmpFileManager.tmp_file
+    with TemporaryDirectory(prefix='sympy_') as tmpdir:
+        with warns(
+            UserWarning,
+            match="The evaluation of the expression is problematic"):
+            i = Integral(log((sin(x)**2 + 1)*sqrt(x**2 + 1)), (x, 0, y))
+            p = plot(i, (y, 1, 5))
+            filename = 'test_advanced_integral.png'
+            p.save(os.path.join(tmpdir, filename))
+            p._backend.close()
+
+
+def test_plot_and_save_5():
+    if not matplotlib:
+        skip("Matplotlib not the default backend")
 
     x = Symbol('x')
     y = Symbol('y')
 
-    p1 = plot(x)
-    p2 = plot_parametric((sin(x), cos(x)), (x, sin(x)), show=False)
-    p3 = plot_parametric(cos(x), sin(x), adaptive=False, nb_of_points=500, show=False)
-    p4 = plot3d_parametric_line(sin(x), cos(x), x, show=False)
-    # symmetric grid
-    p = PlotGrid(2, 2, p1, p2, p3, p4)
-    p.save(tmp_file('%s_grid1' % name))
-    p._backend.close()
-
-    # grid size greater than the number of subplots
-    p = PlotGrid(3, 4, p1, p2, p3, p4)
-    p.save(tmp_file('%s_grid2' % name))
-    p._backend.close()
-
-    p5 = plot(cos(x),(x, -pi, pi), show=False)
-    p5[0].line_color = lambda a: a
-    p6 = plot(Piecewise((1, x > 0), (0, True)), (x, -1, 1), show=False)
-    p7 = plot_contour((x**2 + y**2, (x, -5, 5), (y, -5, 5)), (x**3 + y**3, (x, -3, 3), (y, -3, 3)), show=False)
-    # unsymmetric grid (subplots in one line)
-    p = PlotGrid(1, 3, p5, p6, p7)
-    p.save(tmp_file('%s_grid3' % name))
-    p._backend.close()
-
-
-def test_matplotlib_1():
-
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
-    if matplotlib:
-        try:
-            plot_and_save_1('test')
-        finally:
-            # clean up
-            TmpFileManager.cleanup()
-    else:
-        skip("Matplotlib not the default backend")
-
-def test_matplotlib_2():
-
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
-    if matplotlib:
-        try:
-            plot_and_save_2('test')
-        finally:
-            # clean up
-            TmpFileManager.cleanup()
-    else:
-        skip("Matplotlib not the default backend")
-
-def test_matplotlib_3():
-
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
-    if matplotlib:
-        try:
-            plot_and_save_3('test')
-        finally:
-            # clean up
-            TmpFileManager.cleanup()
-    else:
-        skip("Matplotlib not the default backend")
-
-def test_matplotlib_4():
+    with TemporaryDirectory(prefix='sympy_') as tmpdir:
+        s = Sum(1/x**y, (x, 1, oo))
+        p = plot(s, (y, 2, 10))
+        filename = 'test_advanced_inf_sum.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
-    if matplotlib:
-        try:
-            plot_and_save_4('test')
-        finally:
-            # clean up
-            TmpFileManager.cleanup()
-    else:
-        skip("Matplotlib not the default backend")
+        p = plot(Sum(1/x, (x, 1, y)), (y, 2, 10), show=False)
+        p[0].only_integers = True
+        p[0].steps = True
+        filename = 'test_advanced_fin_sum.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-def test_matplotlib_5():
 
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
-    if matplotlib:
-        try:
-            plot_and_save_5('test')
-        finally:
-            # clean up
-            TmpFileManager.cleanup()
-    else:
+def test_plot_and_save_6():
+    if not matplotlib:
         skip("Matplotlib not the default backend")
 
-def test_matplotlib_6():
+    x = Symbol('x')
 
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
-    if matplotlib:
-        try:
-            plot_and_save_6('test')
-        finally:
-            # clean up
-            TmpFileManager.cleanup()
-    else:
+    with TemporaryDirectory(prefix='sympy_') as tmpdir:
+        filename = 'test.png'
+        ###
+        # Test expressions that can not be translated to np and generate complex
+        # results.
+        ###
+        p = plot(sin(x) + I*cos(x))
+        p.save(os.path.join(tmpdir, filename))
+        p = plot(sqrt(sqrt(-x)))
+        p.save(os.path.join(tmpdir, filename))
+        p = plot(LambertW(x))
+        p.save(os.path.join(tmpdir, filename))
+        p = plot(sqrt(LambertW(x)))
+        p.save(os.path.join(tmpdir, filename))
+
+        #Characteristic function of a StudentT distribution with nu=10
+        x1 = 5 * x**2 * exp_polar(-I*pi)/2
+        m1 = meijerg(((1 / 2,), ()), ((5, 0, 1 / 2), ()), x1)
+        x2 = 5*x**2 * exp_polar(I*pi)/2
+        m2 = meijerg(((1/2,), ()), ((5, 0, 1/2), ()), x2)
+        expr = (m1 + m2) / (48 * pi)
+        p = plot(expr, (x, 1e-6, 1e-2))
+        p.save(os.path.join(tmpdir, filename))
+
+
+def test_plotgrid_and_save():
+    if not matplotlib:
         skip("Matplotlib not the default backend")
 
+    x = Symbol('x')
+    y = Symbol('y')
 
-def test_matplotlib_7():
-
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
-    if matplotlib:
-        try:
-            plotgrid_and_save('test')
-        finally:
-            # clean up
-            TmpFileManager.cleanup()
-    else:
-        skip("Matplotlib not the default backend")
-
+    with TemporaryDirectory(prefix='sympy_') as tmpdir:
+        p1 = plot(x)
+        p2 = plot_parametric((sin(x), cos(x)), (x, sin(x)), show=False)
+        p3 = plot_parametric(
+            cos(x), sin(x), adaptive=False, nb_of_points=500, show=False)
+        p4 = plot3d_parametric_line(sin(x), cos(x), x, show=False)
+        # symmetric grid
+        p = PlotGrid(2, 2, p1, p2, p3, p4)
+        filename = 'test_grid1.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-# Tests for exception handling in experimental_lambdify
-def test_experimental_lambify():
-    x = Symbol('x')
-    f = lambdify([x], Max(x, 5))
-    # XXX should f be tested? If f(2) is attempted, an
-    # error is raised because a complex produced during wrapping of the arg
-    # is being compared with an int.
-    assert Max(2, 5) == 5
-    assert Max(5, 7) == 7
+        # grid size greater than the number of subplots
+        p = PlotGrid(3, 4, p1, p2, p3, p4)
+        filename = 'test_grid2.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
-    x = Symbol('x-3')
-    f = lambdify([x], x + 1)
-    assert f(1) == 2
+        p5 = plot(cos(x),(x, -pi, pi), show=False)
+        p5[0].line_color = lambda a: a
+        p6 = plot(Piecewise((1, x > 0), (0, True)), (x, -1, 1), show=False)
+        p7 = plot_contour(
+            (x**2 + y**2, (x, -5, 5), (y, -5, 5)),
+            (x**3 + y**3, (x, -3, 3), (y, -3, 3)), show=False)
+        # unsymmetric grid (subplots in one line)
+        p = PlotGrid(1, 3, p5, p6, p7)
+        filename = 'test_grid3.png'
+        p.save(os.path.join(tmpdir, filename))
+        p._backend.close()
 
 
 def test_append_issue_7140():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
     if not matplotlib:
         skip("Matplotlib not the default backend")
 
@@ -474,11 +443,8 @@ def test_append_issue_7140():
     with raises(TypeError):
         p1.append(p2._series)
 
-def test_issue_15265():
-    from sympy.core.sympify import sympify
-    from sympy.core.singleton import S
 
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
+def test_issue_15265():
     if not matplotlib:
         skip("Matplotlib not the default backend")
 
@@ -511,27 +477,19 @@ def test_issue_15265():
 
 
 def test_empty_Plot():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
     if not matplotlib:
         skip("Matplotlib not the default backend")
-    from sympy.plotting.plot import Plot
-    p = Plot()
-    # No exception showing an empty plot
-    p.show()
 
-
-def test_empty_plot():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
-    if not matplotlib:
-        skip("Matplotlib not the default backend")
     # No exception showing an empty plot
     plot()
+    p = Plot()
+    p.show()
 
 
 def test_issue_17405():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
     if not matplotlib:
         skip("Matplotlib not the default backend")
+
     x = Symbol('x')
     f = x**0.3 - 10*x**3 + x**2
     p = plot(f, (x, -10, 10), show=False)
@@ -541,9 +499,9 @@ def test_issue_17405():
 
 
 def test_logplot_PR_16796():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
     if not matplotlib:
         skip("Matplotlib not the default backend")
+
     x = Symbol('x')
     p = plot(x, (x, .001, 100), xscale='log', show=False)
     # Random number of segments, probably more than 100, but we want to see
@@ -554,9 +512,9 @@ def test_logplot_PR_16796():
 
 
 def test_issue_16572():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
     if not matplotlib:
         skip("Matplotlib not the default backend")
+
     x = Symbol('x')
     p = plot(LambertW(x), show=False)
     # Random number of segments, probably more than 50, but we want to see
@@ -565,9 +523,9 @@ def test_issue_16572():
 
 
 def test_issue_11865():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
     if not matplotlib:
         skip("Matplotlib not the default backend")
+
     k = Symbol('k', integer=True)
     f = Piecewise((-I*exp(I*pi*k)/k + I*exp(-I*pi*k)/k, Ne(k, 0)), (2*pi, True))
     p = plot(f, show=False)
@@ -578,9 +536,9 @@ def test_issue_11865():
 
 
 def test_issue_11461():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
     if not matplotlib:
         skip("Matplotlib not the default backend")
+
     x = Symbol('x')
     p = plot(real_root((log(x/(x-2))), 3), show=False)
     # Random number of segments, probably more than 100, but we want to see
@@ -588,10 +546,11 @@ def test_issue_11461():
     # and that there are no exceptions.
     assert len(p[0].get_segments()) >= 30
 
+
 def test_issue_11764():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
     if not matplotlib:
         skip("Matplotlib not the default backend")
+
     x = Symbol('x')
     p = plot_parametric(cos(x), sin(x), (x, 0, 2 * pi), aspect_ratio=(1,1), show=False)
     p.aspect_ratio == (1, 1)
@@ -599,13 +558,11 @@ def test_issue_11764():
     # that there are segments generated, as opposed to when the bug was present
     assert len(p[0].get_segments()) >= 30
 
+
 def test_issue_13516():
-    matplotlib = import_module('matplotlib', min_module_version='1.1.0', catch=(RuntimeError,))
-    np = import_module('numpy')
     if not matplotlib:
         skip("Matplotlib not the default backend")
-    if not np:
-        skip("Numpy not found")
+
     x = Symbol('x')
 
     pm = plot(sin(x), backend="matplotlib", show=False)
@@ -623,3 +580,54 @@ def test_issue_13516():
     p = plot(sin(x), show=False)
     assert p.backend == DefaultBackend
     assert len(p[0].get_segments()) >= 30
+
+
+def test_plot_limits():
+    if not matplotlib:
+        skip("Matplotlib not the default backend")
+
+    x = Symbol('x')
+    p = plot(x, x**2, (x, -10, 10))
+    backend = p._backend
+
+    xmin, xmax = backend.ax[0].get_xlim()
+    assert abs(xmin + 10) < 2
+    assert abs(xmax - 10) < 2
+    ymin, ymax = backend.ax[0].get_ylim()
+    assert abs(ymin + 10) < 10
+    assert abs(ymax - 100) < 10
+
+
+def test_plot3d_parametric_line_limits():
+    if not matplotlib:
+        skip("Matplotlib not the default backend")
+
+    x = Symbol('x')
+
+    v1 = (2*cos(x), 2*sin(x), 2*x, (x, -5, 5))
+    v2 = (sin(x), cos(x), x, (x, -5, 5))
+    p = plot3d_parametric_line(v1, v2)
+    backend = p._backend
+
+    xmin, xmax = backend.ax[0].get_xlim()
+    assert abs(xmin + 2) < 1e-2
+    assert abs(xmax - 2) < 1e-2
+    ymin, ymax = backend.ax[0].get_ylim()
+    assert abs(ymin + 2) < 1e-2
+    assert abs(ymax - 2) < 1e-2
+    zmin, zmax = backend.ax[0].get_zlim()
+    assert abs(zmin + 10) < 1e-2
+    assert abs(zmax - 10) < 1e-2
+
+    p = plot3d_parametric_line(v2, v1)
+    backend = p._backend
+
+    xmin, xmax = backend.ax[0].get_xlim()
+    assert abs(xmin + 2) < 1e-2
+    assert abs(xmax - 2) < 1e-2
+    ymin, ymax = backend.ax[0].get_ylim()
+    assert abs(ymin + 2) < 1e-2
+    assert abs(ymax - 2) < 1e-2
+    zmin, zmax = backend.ax[0].get_zlim()
+    assert abs(zmin + 10) < 1e-2
+    assert abs(zmax - 10) < 1e-2