Some years ago I started on something like this; then some years later, I made a package for my differential equations course, in which three plots were aligned to show three different views of equilibria in 1D (autonomous) systems (just as one finds in V. Arnold's Ordinary Differential Equations); some years after that, I started to extend and refactor that package but have not finished. All that is meant to prepare my apology for cutting code out of the package that is complicated and incompletely reworked. It seems to work just fine, which is one reason refactoring it is a low priority.
portrait1D[x'[t] == (1 - x[t]) (x[t] - 3), {x, 0, 4}, {t, 0, 3}]

The code that draws the phase line portrait begins with a (*2*)
in the code dump below and has another comment (*phase line field*)
, which may facilitate searching the long code. Here is the code with appropriate initializations to reproduce the middle plot above:
Block[{
vf = -((-3 + x) (-1 + x)),
fn = x, x1 = 0, x2 = 4,
intervals1D = {{0.`, 1.`}, {1.`, 3.`}, {3.`, 4.`}},
stationaryPts1D = {1, 3},
init1DQ = True,
plotrange, padding},
plotrange = {All, {x1, x2}};
padding = {0.1, (x2 - x1) {0.03, 0.02}};
(*2*)
Show[(*phase line field*)
vectorPlot1D[vf, {fn, x1, x2}, "Vertical",
ImageSize -> {Automatic, 150}],
Graphics[{White, Point[{-0.7, x1}]}] (*extends left margin*),
PlotRange -> plotrange, PlotRangePadding -> padding,
Axes -> {False, True}, AxesLabel -> {None, fn /. x_[0] :> x},
Ticks -> None]
]
There's a polar phase portrait in the package; I could have cut it out, but since it is built on some common utilities, I thought it generous to include it and hope it will be well-received.
polarPortrait1D[
r'[t] == 0.25 r[t] (1 - r[t]) (r[t] - 3),
{r, 0, 4}, {t, 0, 2 Pi}]

(*
* Phase portrait
*)
SetAttributes[odeVFComponents, Listable];
odeVFComponents[de : (Derivative[n_][f_][var_] -> nthDerivative_)] :=
Table[Derivative[i][f][var], {i, n}] /. de;
derivativesToVariables[expr_, fns_, ___] :=
expr /. {Derivative[n0_][f0 : Alternatives @@ Flatten[{fns}]][_] :>
f0[n0], (f1 : Alternatives @@ Flatten[{fns}])[_] :> f1[0]}
doOdeToVF[sols_List] := Flatten@odeVFComponents[sols]
odeToVF[eqn_, fns_, var_] := Module[{solution},
doOdeToVF[solution] /; (
(Head[eqn] === Equal ||
VectorQ[eqn,
MatchQ[#, _Equal] &] || (Message[Solve::naqs, eqn]; False)) &&
Check[
solution = Solve[eqn, Last /@ Sort /@ GatherBy[
Cases[eqn,
Derivative[_][Alternatives @@ Flatten[{fns}]][var],
Infinity]
, # /. Derivative[_] -> Identity &]
]; True, False, {Solve::ivar}]
)
];
(* need to convert ODE, List VFs to scalar VFs
* to map var x[t] to x if necessary
* to check vf is 1d func of x only
* to initialize singularities and intervals
* to add stationary point to de object
*)
ClearAll[init1D, vectorPlot1D, portrait1D, polarPortrait1D];
SetAttributes[init1D, HoldAll];
(*
* This is an important hook to call init1D on any 1D system plotter
* init1D[] parses the ode/vector field and sets up stationary \
solutions
* and the intervals between them
*)
$plotters1D = vectorPlot1D | portrait1D | polarPortrait1D;
(call : #[args___] /; ! TrueQ[init1DQ] :=
init1D[call]) & /@ $plotters1D;
(* * * * *)
zeros[vf_, {var_, a_, b_}] :=
var /. Solve[vf == 0 && a <= var <= b, var, Reals] /. var -> {};
parseODE1D[ode_Equal, {x_[t_] | x_, x1_, x2_}, T : {t_, t1_, t2_},
rest___] :=
{odeToVF[ode, {x}, t] /. x[t] -> x, {x, x1, x2}, T,
rest}
parseODE1D[ode_Equal, {x_[t_], x1_, x2_}, rest___] /; !
NumericQ[t] :=
{odeToVF[ode, {x}, t] /. x[t] -> x, {x, x1, x2},
rest};
parseODE1D[de : diffEq[asc_], rest___] := With[{
x = getKey[de, "dependentVars", {}] /. {t1_} :> t1,
t = getKey[de, "independentVars", {}] /. {t1_} :> t1},
{odeToVF[getKey[de, "de", {}], {x}, t] /. x[t] -> x, rest} (**
TBD: Check rest **)
];
parseODE1D[___] := TBD["input error (diagnosis TBD)"]; (** TBD ??? **)
init1D[h_[ode : _Equal | _diffEq, args___]] :=
init1D[h[##]] & @@ parseODE1D[ode, args];
init1D[h_[vf_, X : {x_, a_, b_}, args___]] :=
Block[{init1DQ = True, stationaryPts1D, intervals1D}
, stationaryPts1D = zeros[#, X]
; intervals1D =
Partition[Union[N[Join[{a, b}, stationaryPts1D]]], 2, 1]
; h[#, X, args]
] & /@ Flatten@vf /. {res_} :> res; (** Flatten@vf ??? **)
vectorPlot1D // Options = Join[{
VectorPoints -> Automatic,
PlotStyle -> Automatic,
"VectorColors" -> {ColorData[1][1], ColorData[1][2]}},
Options[Graphics]];
(* new diffEqClass` interfaces *)
vectorPlot1D[vf_, X : {_, _, _}, T : {_, _, _}, rest___] :=
vectorPlot1D[vf, X, rest];
vectorPlot1D[{{vf_}}, {var_, a_, b_}, (* Needed ??? *)
Optional[orientation : "Horizontal" | "Vertical", "Vertical"],
opts : OptionsPattern[]] /; True :=
vectorPlot1D[vf, {var, a, b}, orientation, opts];
vectorPlot1D[vfs : {Repeated[{_}, {2, Infinity}]},
Optional[orientation : "Horizontal" | "Vertical",
"Vertical"], {var_, a_, b_}, opts : OptionsPattern[]] :=
vectorPlot1D[#, {var, a, b}, orientation, opts] & /@ Flatten@vfs;
(*vectorPlot1D[bad_,{var_,a_,b_},opts:OptionsPattern[]]:=(
Message[diffEq::CantHappen,"Plotting",Row[{"the 1D vector field is \
",bad}']];
Throw[$Failed]);*)
(* old math212` function *)
vectorPlot1D[vf_, {var_, a_, b_},
Optional[orientation : "Horizontal" | "Vertical", "Vertical"],
opts : OptionsPattern[]] :=
Module[{style, colorFn, coord}
, If[orientation === "Horizontal"
, coord[y_][x_] := {x, y}
, coord[x_][y_] := {x, y}]
; style = OptionValue[PlotStyle] /.
Automatic -> Directive[Thickness[Medium], ColorData[1][1]]
; colorFn[{y1_, y2_}] :=
If[(vf /. var -> Mean[{y1, y2}]) > 0
, Sequence @@ OptionValue["VectorColors"]]
; Graphics[{
Thickness[Medium], ColorData[1][1], style
, {colorFn[#1], Line[Thread[coord[0][#]]]} & /@ intervals1D
, {colorFn[#1], Arrowheads[Medium]
, Arrow[If[(vf /. var -> Mean[#1]) > 0
, Identity
, Reverse
][Thread[coord[0.08 (b - a)][With[{scale = 0.9},
{{scale, 1 - scale}, {1 - scale, scale}} . #1]]]
]
]
} & /@ SortBy[N[intervals1D], Differences]
, EdgeForm[style]
, {Black, Text[Round[N[#1], 0.01], coord[-0.1][#], coord[1][0]]
, White, Disk[coord[0][#], Offset[{2, 2}]]
} & /@ stationaryPts1D
}
, FilterRules[{opts}, Options[Graphics]]
]
];
(* is fn/.x_[0]\[RuleDelayed]x still necessary? useful? *)
portrait1D[vf_, {fn_, x1_, x2_}, {time_, t1_, t2_}, opts___?OptionQ] :=
Module[{plotrange = {All, {x1, x2}},
padding = {0.1, (x2 - x1) {0.03, 0.02}}}
, Grid[{
{(*1*)
ParametricPlot[ (* velocity plot *)
{vf, fn}, {fn, x1, x2}
, Mesh -> {{0}}, MeshFunctions -> {#1 &}, Ticks -> None,
ImageSize -> {Automatic, 150}, PlotRange -> plotrange,
PlotRangePadding -> padding , AspectRatio -> 2,
AxesLabel -> {HoldForm[Global`v], fn /. x_[0] :> x},
Epilog -> {White, EdgeForm[Black],
Disk[{0, #}, Offset[{2, 2}]] & /@ stationaryPts1D}
]
,(*2*)
Show[ (* phase line field *)
vectorPlot1D[vf, {fn, x1, x2}, ImageSize -> {Automatic, 150}]
, Graphics[{White, Point[{-0.7, x1}]}] (*
extends left margin *)
, PlotRange -> plotrange,
PlotRangePadding -> padding, Axes -> {False, True},
AxesLabel -> {None, fn /. x_[0] :> x}, Ticks -> None
]
,(*3*)
Module[{eqns, x0} (* typical solutions *)
,
eqns = {fn'[time] == (vf /. fn -> fn[time]), fn[0] == x0}
; With[{psol = ParametricNDSolveValue[eqns
, fn, {time, t1, t2}, {x0}
,
"ExtrapolationHandler" -> {Indeterminate &,
"WarningMessage" -> False}
]}
, Plot[Evaluate[Quiet[
(psol[#1][time] &) /@
Join[stationaryPts1D
,
If[(vf /. fn -> Mean[#1]) >
0, {0.9, 0.1} . #1, {0.1, 0.9} . #1] & /@
intervals1D
]]]
, {time, t1, t2}
, PlotRange -> plotrange,
PlotRangePadding -> padding , ImageSize -> {Automatic, 150},
Ticks -> None, ExclusionsStyle -> None,
AxesLabel -> {time, fn /. x_[0] :> x},
Epilog -> {White, EdgeForm[Black],
Disk[{0, #}, Offset[{2, 2}]] & /@ stationaryPts1D}
]
]
]
},
(Row[#1, BaseStyle -> "Label"] &) /@ {{fn /. x_[0] :> x, " vs. ",
HoldForm[Global`v]}, {"vector field"}, {fn /. x_[0] :> x,
" vs. " , time}}
}]
];
(*
* Polar direction field portrait
* alternative APIs TBD ???
*)
polarPortrait1D[vf_, {fn_, r1_, r2_}, {time_, t1_, t2_},
opts : OptionsPattern@PolarPlot] :=
Block[{plotrange = r2, padding = {0.1, (r2 - r1) {0.03, 0.02}}}
, Module[{eqns, x0}
, eqns = {
fn'[time] == (vf /. fn -> fn[time]),
fn[t1] == x0, (* IC *)
Table[
With[{r0 = r0},
WhenEvent[Abs[fn[time] - r0] < 1*^-3, "StopIntegration"]],
{r0, stationaryPts1D}]}
; With[
{psol = ParametricNDSolveValue[eqns, fn, {time, t1, t2}, {x0},
"ExtrapolationHandler" -> {Indeterminate &,
"WarningMessage" -> False}]}
, PolarPlot[
Evaluate@Flatten[List[
ConditionalExpression[#, t1 <= time <= t1 + 2 Pi] & /@
DeleteDuplicates@stationaryPts1D,
Quiet[(psol[#1]@(time) &) /@ (If[(vf /. fn -> Mean[#1]) >
0, {0.9, 0.1} . #1, {0.1, 0.9} . #1] &) /@
intervals1D]]],
{time, t1, t2},
opts,
PlotRange -> plotrange, PlotRangePadding -> padding,
Ticks -> None, ExclusionsStyle -> None,
PlotLabel -> First@eqns] /.
l_Line :> {Arrowheads[
PadRight[{0.}, 1 + Ceiling[2/r2*ArcLength@l], 0.04]],
Arrow @@ l}
]
]
];