diff --git a/doc/generic/pgf/CHANGELOG.md b/doc/generic/pgf/CHANGELOG.md
index a9d8f9ade..37c72b9be 100644
--- a/doc/generic/pgf/CHANGELOG.md
+++ b/doc/generic/pgf/CHANGELOG.md
@@ -11,6 +11,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 ### Added
 
 - Add `RGB` and `gray` color model support for ConTeXt #1130
+- Add a new monotonic interpolation plot handler #1358
 
 ### Fixed
 
diff --git a/doc/generic/pgf/pgfmanual-en-library-plot-handlers.tex b/doc/generic/pgf/pgfmanual-en-library-plot-handlers.tex
index 689e3f0c1..e8fe5268d 100644
--- a/doc/generic/pgf/pgfmanual-en-library-plot-handlers.tex
+++ b/doc/generic/pgf/pgfmanual-en-library-plot-handlers.tex
@@ -97,6 +97,38 @@ \subsection{Curve Plot Handlers}
 \end{command}
 
 
+\begin{command}{\pgfplothandlermonotone}
+    This handler will issue a |\pgfpathcurveto| command for each point of the
+    plot, \emph{except} possibly for the first. As for the line-to handler,
+    what happens with the first point can be specified using
+    |\pgfsetmovetofirstplotpoint| or |\pgfsetlinetofirstplotpoint|.
+
+    Obviously, the |\pgfpathcurveto| command needs, in addition to the points
+    on the path, some control points. These are computed using The
+    Fritsch-Carlson monotone cubic interpolation algorithm which roughly does
+    the following:
+
+    \begin{enumerate}
+        \item Initialize tangent slopes at each input point as the
+            average of the slopes of the neighboring secant lines.
+        \item If the secants have opposite signes (because this is a local
+            extremum), set the tangent slope to~$0$ to get an horizontal
+            tangent and prevent overshoot. Same if one of the neighboring
+            secants is horizontal.
+        \item Adjust the tangent slopes to ensure strict monotonicity.
+        \item Place the control points on the computed tangent lines,
+            at one third of the $x$~distance of a secant line: if $a$, $b$
+            and~$c$ are three consecutive input points and $m$~is the
+            slope computed at~$B$, then the control points around~$b$ will
+            be $u$~and~$v$ where $x_u = x_b - \frac{1}{3}(x_b-x_a)$,
+            $y_u = y_b - \frac{1}{3} m (x_b-x_a)$,
+            $x_v = x_b + \frac{1}{3}(x_c-x_b)$ and
+            $y_v = y_b + \frac{1}{3} m (x_c-x_b)$.
+    \end{enumerate}
+\end{command}
+
+
+
 \subsection{Constant Plot Handlers}
 
 There are several plot handlers which produce piecewise constant interpolations
diff --git a/tex/generic/pgf/frontendlayer/tikz/libraries/datavisualization/tikzlibrarydatavisualization.code.tex b/tex/generic/pgf/frontendlayer/tikz/libraries/datavisualization/tikzlibrarydatavisualization.code.tex
index 04dae1e6b..a39ae8fc9 100644
--- a/tex/generic/pgf/frontendlayer/tikz/libraries/datavisualization/tikzlibrarydatavisualization.code.tex
+++ b/tex/generic/pgf/frontendlayer/tikz/libraries/datavisualization/tikzlibrarydatavisualization.code.tex
@@ -1538,6 +1538,7 @@
   smooth cycle/.style={@set={\pgfplothandlerclosedcurve}{default label in legend closed path}},
   straight line/.style={@set={\pgfplothandlerlineto}{default label in legend path}},
   straight cycle/.style={@set={\pgfplothandlerpolygon}{default label in legend closed path}},
+  smooth monotone line/.style={@set={\pgfplothandlermonotone}{default label in legend path}},
   polygon/.style={straight cycle},% alias
   gap line/.style={@set={\pgfplothandlergaplineto}{default label in legend path}},
   gap cycle/.style={@set={\pgfplothandlergapcycle}{gap circular label in legend line}}
diff --git a/tex/generic/pgf/frontendlayer/tikz/tikz.code.tex b/tex/generic/pgf/frontendlayer/tikz/tikz.code.tex
index d9b1743b9..957dfc23c 100644
--- a/tex/generic/pgf/frontendlayer/tikz/tikz.code.tex
+++ b/tex/generic/pgf/frontendlayer/tikz/tikz.code.tex
@@ -1232,6 +1232,7 @@
 % Plot options
 \tikzoption{smooth}[]{\let\tikz@plot@handler=\pgfplothandlercurveto}%
 \tikzoption{smooth cycle}[]{\let\tikz@plot@handler=\pgfplothandlerclosedcurve}%
+\tikzoption{smooth monotone}[]{\let\tikz@plot@handler=\pgfplothandlermonotone}%
 \tikzoption{sharp plot}[]{\let\tikz@plot@handler\pgfplothandlerlineto}%
 \tikzoption{sharp cycle}[]{\let\tikz@plot@handler\pgfplothandlerpolygon}%
 
diff --git a/tex/generic/pgf/libraries/pgflibraryplothandlers.code.tex b/tex/generic/pgf/libraries/pgflibraryplothandlers.code.tex
index d0b9a388e..dde258b49 100644
--- a/tex/generic/pgf/libraries/pgflibraryplothandlers.code.tex
+++ b/tex/generic/pgf/libraries/pgflibraryplothandlers.code.tex
@@ -269,6 +269,135 @@
 }%
 
 
+% This handler converts each plot stream command into a curveto
+% command, except for the first, which is converted to the previously
+% specified action. The curveto control points are computed and adjusted
+% so that the resulting curve respects the monotonicity of the points.
+% If a point is a local extremum of the point set, it will also be
+% a local extremum of the generated curve.
+% This handler treats the x and y direction differently, and expects
+% x values to be monotone.
+% See https://en.wikipedia.org/wiki/Monotone_cubic_interpolation for an
+% explanation of the algorithm. The last part uses \pgfpathcurveto
+% instead of using Hermitte base functions.
+\pgfdeclareplothandler{\pgfplothandlermonotone}{}{%
+    point macro=\pgf@plot@monotone@dopointone,
+    jump macro=\pgf@plot@monotone@dojump,
+    end macro=\pgf@plot@monotone@doend
+}%
+
+\def\pgf@plot@monotone@dojump{%
+    \pgf@plot@monotone@doend
+    \global\pgf@plot@startedfalse
+    \global\let\pgf@plotstreampoint\pgf@plot@monotone@dopointone
+}%
+
+\def\pgf@plot@monotone@dopointone#1{%
+    \pgf@process{#1}%
+    \pgf@plot@first@action{\pgfqpoint{\pgf@x}{\pgf@y}}%
+    \xdef\pgf@plot@monotone@pointone{\noexpand\pgfqpoint{\the\pgf@x}{\the\pgf@y}}%
+    \global\let\pgf@plotstreampoint=\pgf@plot@monotone@dopointtwo
+}%
+
+\def\pgf@plot@monotone@dopointtwo#1{%
+    \pgf@process{#1}%
+    \xdef\pgf@plot@monotone@pointtwo{\noexpand\pgfqpoint{\the\pgf@x}{\the\pgf@y}}%
+    % compute vector:
+    \pgf@xa=\pgf@x%
+    \pgf@ya=\pgf@y%
+    \pgf@process{\pgf@plot@monotone@pointone}%
+    \advance\pgf@xa by-\pgf@x
+    \advance\pgf@ya by-\pgf@y
+    % save left delta
+    \xdef\pgf@plot@monotone@deltaleft{\pgf@sys@tonumber\pgf@xa}%
+    % compute slope
+    \pgfmathdivide@{\the\pgf@ya}{\the\pgf@xa}%
+    \global\let\pgf@plot@monotone@rateleft\pgfmathresult
+    % init tangent one
+    \global\let\pgf@plot@monotone@slopeone\pgf@plot@monotone@rateleft
+    % prepare for next step
+    \global\let\pgf@plotstreampoint=\pgf@plot@monotone@dopointthree
+    \global\pgf@plot@startedtrue%
+}%
+
+\def\pgf@plot@monotone@dopointthree#1{%
+    \pgf@process{#1}%
+    \xdef\pgf@plot@monotone@pointthree{\noexpand\pgfqpoint{\the\pgf@x}{\the\pgf@y}}%
+    % compute vector:
+    \pgf@xa=\pgf@x%
+    \pgf@ya=\pgf@y%
+    \pgf@process{\pgf@plot@monotone@pointtwo}%
+    \advance\pgf@xa by-\pgf@x
+    \advance\pgf@ya by-\pgf@y
+    % save right delta
+    \xdef\pgf@plot@monotone@deltaright{\pgf@sys@tonumber\pgf@xa}%
+    % compute slope
+    \pgfmathdivide@{\the\pgf@ya}{\the\pgf@xa}%
+    \global\let\pgf@plot@monotone@rateright\pgfmathresult
+    % compute tangent slope
+    \pgfmathadd@{\pgf@plot@monotone@rateleft}{\pgf@plot@monotone@rateright}%
+    \pgf@xa=\pgfmathresult pt\relax
+    \pgf@xa=0.5\pgf@xa
+    % Fix slopes & draw curve
+    \pgf@plot@monotonecurveto@fixdraw
+    % Prepare next:
+    \global\let\pgf@plot@monotone@slopeone\pgf@plot@monotone@slopetwo
+    \global\let\pgf@plot@monotone@pointone\pgf@plot@monotone@pointtwo
+    \global\let\pgf@plot@monotone@pointtwo\pgf@plot@monotone@pointthree
+    \global\let\pgf@plot@monotone@deltaleft\pgf@plot@monotone@deltaright
+    \global\let\pgf@plot@monotone@rateleft\pgf@plot@monotone@rateright
+}%
+
+\def\pgf@plot@monotonecurveto@fixdraw{%
+    % fix extremums
+    % \pgf@xa is slope two
+    \pgf@ya=\pgf@plot@monotone@rateleft pt\relax
+    \pgf@yb=\pgf@plot@monotone@rateright\pgf@ya
+    \ifdim\pgf@yb<0.0001pt\relax
+        \pgf@xa=0pt\relax
+    \else
+    % fix motonicity
+        \pgf@ya=3\pgf@ya
+        \ifdim\pgf@xa>0pt\relax
+            \ifdim\pgf@xa>\pgf@ya \pgf@xa=\pgf@ya\fi
+        \else
+            \ifdim\pgf@xa<\pgf@ya \pgf@xa=\pgf@ya\fi
+        \fi
+        \pgf@ya=\pgf@plot@monotone@rateright pt\relax
+        \pgf@ya=3\pgf@ya
+        \ifdim\pgf@xa>0pt\relax
+            \ifdim\pgf@xa>\pgf@ya \pgf@xa=\pgf@ya\fi
+        \else
+            \ifdim\pgf@xa<\pgf@ya \pgf@xa=\pgf@ya\fi
+        \fi
+    \fi
+    \xdef\pgf@plot@monotone@slopetwo{\pgf@sys@tonumber\pgf@xa}%
+    % compute points for left bezier
+    % (\pgf@x,\pgf@y) is point two
+    \pgf@ya=\pgf@plot@monotone@deltaleft pt\relax
+    \pgf@ya=0.333333\pgf@ya
+    \advance \pgf@x by -\pgf@ya
+    \advance \pgf@y by -\pgf@plot@monotone@slopetwo\pgf@ya
+    \xdef\pgf@plot@monotone@supporttwo{\noexpand\pgfqpoint{\the\pgf@x}{\the\pgf@y}}%
+    \pgf@process{\pgf@plot@monotone@pointone}%
+    \advance \pgf@x by \pgf@ya
+    \advance \pgf@y by \pgf@plot@monotone@slopeone\pgf@ya
+    \xdef\pgf@plot@monotone@supportone{\noexpand\pgfqpoint{\the\pgf@x}{\the\pgf@y}}%
+    % draw left curve:
+    \pgfpathcurveto{\pgf@plot@monotone@supportone}{\pgf@plot@monotone@supporttwo}{\pgf@plot@monotone@pointtwo}%
+}%
+
+\def\pgf@plot@monotone@doend{%
+    \ifpgf@plot@started
+        % fixdraw needs (\pgf@x,\pgf@y) as point two
+        % and \pgf@xa as slope two
+        \pgf@process{\pgf@plot@monotone@pointtwo}%
+        \pgf@xa=\pgf@plot@monotone@rateleft pt\relax
+        \pgf@plot@monotonecurveto@fixdraw
+    \fi
+}%
+
+