How to plot logistic regression decision boundary?
Your decision boundary is a surface in 3D as your points are in 2D.

With Wolfram Language

Create the data sets.

mqtrue = 5;
cqtrue = 30;
With[{x = Subdivide[0, 3, 50]},
  dat1 = Transpose@{x, mqtrue x + 5 RandomReal[1, Length@x]};
  ];
With[{x = Subdivide[7, 10, 50]},
  dat2 = Transpose@{x, mqtrue x + cqtrue + 5 RandomReal[1, Length@x]};
  ];


View in 2D (ListPlot) and the 3D (ListPointPlot3D) regression space.

ListPlot[{dat1, dat2}, PlotMarkers -&gt; "OpenMarkers", PlotTheme -&gt; "Detailed"]



  


I Append the response variable to the data. 

datPlot =
 ListPointPlot3D[
  MapThread[Append, {#, Boole@Thread[#[[All, 2]] &gt; 40]}] &amp; /@ {dat1, dat2}
  ]



  


Perform a Logistic regression (LogitModelFit). You could use GeneralizedLinearModelFit with ExponentialFamily set to "Binomial" as well.

With[{dat = Join[dat1, dat2]},
 model =
  LogitModelFit[
   MapThread[Append, {dat, Boole@Thread[dat[[All, 2]] &gt; 40]}],
   {x, y}, {x, y}]
 ]



  


From the FittedModel "Properties" we need "Function".

model["Properties"]



  {AdjustedLikelihoodRatioIndex,    DevianceTableDeviances, ParameterConfidenceIntervalTableEntries,
  AIC,    DevianceTableEntries,   ParameterConfidenceRegion,
  AnscombeResiduals,  DevianceTableResidualDegreesOfFreedom,  ParameterErrors,
  BasisFunctions, DevianceTableResidualDeviances, ParameterPValues,
  BestFit,    EfronPseudoRSquared,    ParameterTable,
  BestFitParameters,  EstimatedDispersion,    ParameterTableEntries,
  BIC,    FitResiduals,   ParameterZStatistics,
  CookDistances,  Function,   PearsonChiSquare,
  CorrelationMatrix,  HatDiagonal,    PearsonResiduals,
  CovarianceMatrix,   LikelihoodRatioIndex,   PredictedResponse,
  CoxSnellPseudoRSquared, LikelihoodRatioStatistic,   Properties,
  CraggUhlerPseudoRSquared,   LikelihoodResiduals,    ResidualDeviance,
  Data,   LinearPredictor,    ResidualDegreesOfFreedom,
  DesignMatrix,   LogLikelihood,  Response,
  DevianceResiduals,  NullDeviance,   StandardizedDevianceResiduals,
  Deviances,  NullDegreesOfFreedom,   StandardizedPearsonResiduals,
  DevianceTable,  ParameterConfidenceIntervals,   WorkingResiduals,
  DevianceTableDegreesOfFreedom,  ParameterConfidenceIntervalTable}   


model["Function"]



  


Use this for prediction 

model["Function"][8, 54]



0.0196842



and plot the decision boundary surface in 3D along with the data (datPlot) using Show and Plot3D

modelPlot =
 Show[
  datPlot,
  Plot3D[
   model["Function"][x, y],
   Evaluate[
    Sequence @@ 
     MapThread[Prepend, {MinMax /@ Transpose@Join[dat1, dat2], {x, y}}]],
   Mesh -&gt; None,
   PlotStyle -&gt; Opacity[.25, Green],
   PlotPoints -&gt; 30
   ]
  ]




With ParametricPlot3D and Manipulate you can examine decision boundary curves for values of the variables. For example, keeping x fixed and letting y vary or vice versa.

Manipulate[
 Show[
  modelPlot,
  ParametricPlot3D[
   {x, u, model["Function"][x, u]}, {u, 0, 80}, PlotStyle -&gt; Orange],
  ParametricPlot3D[
   {u, y, model["Function"][u, y]}, {u, 0, 10}, PlotStyle -&gt; Purple],
  PlotLabel -&gt; 
   StringTemplate["model[`1`, `2`] = `3`"] @@ {x, y, model["Function"][x, y]}
  ],
 {{x, 6, Style["x", Orange, Bold]}, 0, 10, Appearance -&gt; "Labeled"},
 {{y, 40, Style["y", Purple, Bold]}, 0, 80, Appearance -&gt; "Labeled"}
 ]




Update

You can also plot contours of the probability in 2D.

plot = ListPlot[{dat1, dat2}, PlotMarkers -&gt; "OpenMarkers", PlotTheme -&gt; "Detailed"];

Manipulate[
 db = y /. First@Quiet@Solve[model["Function"][x, y] == p, y];
 Show[
  plot,
  Plot[db, {x, 0, 10}, PlotStyle -&gt; Red]
  ],
 {{p, .5}, 0, 1, Appearance -&gt; "Labeled"}
 ]




Hope this helps.