You want an example for solving least absolute deviation by linear programming. I will show you an simple implementation in R. Quantile regression is a generalization of least absolute deviation, which is the case of the quantile 0.5, so I will show a solution for quantile regression. Then you can check the results with the R quantreg package:

    rq_LP  <-  function(x, Y, r=0.5, intercept=TRUE) {
        require("lpSolve")
        if (intercept) X  <-  cbind(1, x) else X <-  cbind(x)
        N   <-  length(Y)
        n  <-  nrow(X)
        stopifnot(n == N)
        p  <-  ncol(X)
        c  <-  c(rep(r, n), rep(1-r, n), rep(0, 2*p))  
                 # cost coefficient vector
        A  <- cbind(diag(n), -diag(n), X, -X)
        res  <-  lp("min", c, A, "=", Y, compute.sens=1)
            ### Desempaquetar los coefs:
        sol <- res$solution
        coef1  <-  sol[(2*n+1):(2*n+2*p)]
        coef <- numeric(length=p)
        for (i in seq(along=coef)) {
             coef[i] <- (if(coef1[i]<=0)-1 else +1) *  
                  max(coef1[i], coef1[i+p])
        }
        return(coef)
        }

Then we use it in a simple example:

    library(robustbase)
    data(starsCYG)
    Y  <- starsCYG[, 2]
    x  <- starsCYG[, 1]
    rq_LP(x, Y)
    [1]  8.1492045 -0.6931818

then you yourself can do the check with quantreg.

I've seen a comment saying that I can't treat the absolute value the way I do, but the logic is as follows.

We have e_i = | y_i - a - b x_i |

Since we are minimising the sum of e_i, at the optimal point, this is equivalent to e_i >= | y_i - a -b x_i |. From there, it is separated into positive and negative cases.

Let's formulate the Sum of Absolute Residuals / Least Absolute Deviation (LAD) / $ {L}_{1} $ Regression / Robust Regression:

$$ \arg \min_{x} {\left\| A x - b \right\|}_{1} $$

There are few options to solve this.

Iteratively Reweighted Least Squares (IRLS)

The main trick in the IRLS approach is that:

$$ {\left\| \boldsymbol{x} \right\|}_{p}^{p} = \sum_{i} {\left| {x}_{i} \right|}^{p - 2} {\left| {x}_{i} \right|}^{2} = \sum_{i} {w}_{i}^{2} {\left| {x}_{i} \right|}^{2} $$

On the other hand, it is known (Weighted Least Squares):

$$ \sum_{i = 1}^{m} {w}_{i} {\left| {A}_{i, :} x - {b}_{i} \right|}^{2} = {\left\| {W}^{\frac{1}{2}} \left( A x - b \right) \right\|}_{2}^{2} $$

The solution is given by:

$$ \arg \min_{x} {\left\| {W}^{\frac{1}{2}} \left( A x - b \right) \right\|}_{2}^{2} = {\left( {A}^{T} W A \right)}^{-1} {A}^{T} W b $$

In the case above, where $ p = 1 $, one could set $ {w}_{i} = {\left| {A}_{i, :} x - {b}_{i} \right|}^{-1} $.
Since the solution depends on $ x $ one could set iterative procedure:

$$ {x}^{k + 1} = {\left( {A}^{T} {W}^{k} A \right)}^{-1} {A}^{T} {W}^{k} b $$

Where $ {W}_{ii}^{k} = {w}_{i}^{k} = {\left| {A}_{i, :} {x}^{k} - {b}_{i} \right|}^{p - 2} = {\left| {A}_{i, :} {x}^{k} - {b}_{i} \right|}^{-1} $.

Sub Gradient Method

The Sub Gradient is given by:

$$ \frac{\partial {\left\| A x - b \right\|}_{1}}{\partial x} = {A}^{T} \operatorname{sgn} \left( A x - b \right) $$

So the solution is given by iterating:

$$ {x}^{k + 1} = {x}^{k} - \alpha \left[ {A}^{T} \operatorname{sgn} \left( A x - b \right) \right] $$

Linear Programming Conversion

This approach convert the problem into Linear Programming problem.

By defining $ {t}_{i} = \left| {A}_{i, :} x - {b}_{i} \right| $ the optimization problem becomes:

$$\begin{align*} \text{minimize} \quad & \boldsymbol{1}^{T} \boldsymbol{t} \\ \text{subject to} \quad & -\boldsymbol{t} \preceq A x -b \preceq \boldsymbol{t} \\ & \boldsymbol{0} \preceq \boldsymbol{t} \end{align*}$$

Which is a Linear Programming form of the problem.

MATLAB Code

I implemented and validated the approaches by comparing it to CVX solution. The MATLAB code is available at my Mathematics StackExchange Question 2603548 Repository.

Remark: Pay attention that the MATLAB code is a vanilla implementation. Real world implementation should take into account numerical issues (Like the inverse of a small number).