For Wolfram Language you may use ImplicitRegion, ContourPlot, and RegionPlot.
For a system
system =
{
y >= 2
, 2 y <= 25 - x
, 4 y >= 2 x - 8
, y <= 2 x - 5
};
The ImplicitRegion and its RegionBounds (useful for plotting) can be obtained with
fregion = ImplicitRegion[system, {x, y}];
rbounds = RegionBounds@fregion;
plotbounds = Transpose[{.9, 1.1} Transpose@rbounds];
fregion can be used with all Region Properties and Measures if you need addtional information on the feasible region.
The inequalities of system can be ContourPloted by Applying Equal. Evaluate is used to resolve the expressions before the passing them to ContourPlot; only needed because converting the inequalities and determining the bounds on the fly.
frplot =
Show[
ContourPlot[
Evaluate[Equal @@@ system]
, Evaluate[Sequence @@ MapThread[Prepend, {plotbounds, {x, y}}]]
, PlotLegends -> system
, AspectRatio -> Automatic
]
, RegionPlot[fregion, BoundaryStyle -> None]
]

We can check a couple of objection function solutions.
objectives = {x y, x - 2 y};
sols = Maximize[{#, system}, {x, y}] & /@ objectives
{{625/8, {x->25/2, y->25/4}}, {4, {x->8, y->2}}}
ListPlot with Callout labels.
splot =
ListPlot[
MapIndexed[
Callout[{x, y} /. Last@#, "Max" <> ToString@objectives[[#2]]] &,
sols]
]

Then Show the plots together.
Show[frplot, splot]

Hope this helps.