This is the Affine Scaling Algorithm outlined by:

Linear and nonlinear programming with Maple: an interactive, applications ...
Paul E. Fishback,Paul F. Fishback

http://books.google.se/books?id=qwRLD9oXvQAC&pg=PA68&lpg=PA68&dq=%22Linear+and+Nonlinear+Programming+with+Maple%22+Section+2.6:+Interior+Point+Algorithm+I&source=bl&ots=gvpyCxR-eA&sig=zN00EKB-syghqS0zK2xuEa9rdXk&hl=sv&ei=nLGuToDuBOT04QSy1-n4Dg&sa=X&oi=book_result&ct=result&resnum=3&ved=0CCsQ6AEwAg#v=onepage&q=Fuel%20LP%20Pro&f=false


restart:
with(plots):
with(Optimization):

Ob := 4*x[1]+3*x[2]:
con1 := x[1] <= 8:
con2 := 2*x[1]+2*x[2] <= 28:
con3 := 3*x[1]+2*x[2] <= 32:
con4 := x[1] >= 0:
con5 := x[2] >= 0:

con := [con1, con2, con3, con4, con5]:

evalf(simplex[maximize](Ob, con), 3);
Maximize(Ob, con, iterationlimit = 1000);
constraintsA := {con1, con2, con3, con4, con5};

inequal(constraintsA, x[1] = 0 .. 15, x[2] = 0 .. 20, labels = ["x[1]", "x[2]"], labelfont = [TIMES, ROMAN, 14]);


                                         {x[1] = 4., x[2] = 10.}
             [46., [x[1] = 4.00000000000000178, x[2] = 9.9999999999999982]]
  constraintsA := {0 <= x[1], 0 <= x[2], x[1] <= 8, 2 x[1] + 2 x[2] <= 28, 3 x[1] + 2 x[2] <= 32}





restart:
with(plots):
with(LinearAlgebra):

c := Vector[row]([4, 3, 0, 0, 0]):
X := Vector[row]([x[1], x[2], x[3], x[4], x[5]]):
c.X:
A := Matrix(3, 5, [1, 0, 1, 0, 0, 2, 2, 0, 1, 0, 3, 2, 0, 0, 1]):
b := `<,>`(8, 28, 32):
N := 15:
lambda := .75:
x := array(0 .. N):
x[0] := `<,>`(3, 4, 5, 14, 15):

for i from 0 to N-1 do
d := DiagonalMatrix(convert(x[i], list));
`~A` := A.d;
`~c` := c.d;
P := IdentityMatrix(5)-Transpose(`~A`).MatrixInverse(`~A`.Transpose(`~A`)).`~A`;
cp := P.Transpose(`~c`);
alpha := 1/abs(min(seq(cp[j], j = 1 .. 5)));
x[i+1] := d.(`<,>`(1, 1, 1, 1, 1)+alpha*lambda.cp)
end do:

interface(rtablesize = 20):
Matrix([convert([0, init, seq(it[i], i = 1 .. N-1)], Vector), Transpose(Matrix([convert([x1, x2, slack1, slack2, slack3, OB], Vector), seq(`<,>`(x[i], 4*x[i][1]+3*x[i][2]), i = 0 .. N-1)]))]);

constraintsA := {xx[1] >= 0, xx[2] >= 0, xx[1] <= 8, 2*xx[1]+2*xx[2] <= 28, 3*xx[1]+2*xx[2] <= 32}:

Ap := Array(1 .. 1, 1 .. 2):
Ap[1, 1] := display({inequal(constraintsA, xx[1] = 0 .. 15, xx[2] = 0 .. 20, labels = ["x[1]", "x[2]"], labelfont = [TIMES, ROMAN, 14]), pointplot([seq([x[i][1], x[i][2]], i = 0 .. N-1)], color = red, thickness = 3, connect = true), pointplot([seq([x[i][1], x[i][2]], i = 0 .. N-1)], color = red, thickness = 3, symbol = solidcircle, symbolsize = 14)}, title = "Constraint", titlefont = [TIMES, ROMAN, 16]):
Ap[1, 2] := display({pointplot([seq([i, 4*x[i][1]+3*x[i][2]], i = 0 .. N-1)], connect = true), pointplot([seq([i, 4*x[i][1]+3*x[i][2]], i = 0 .. N-1)], color = black, thickness = 3, symbol = solidcircle, symbolsize = 14)}, title = "Objective Function", titlefont = [TIMES, ROMAN, 16]):
display(Ap)