Question: LP and Interior Point Method

I have been watching a demonstration on youtube about LP and Interior Point Methods:

http://www.youtube.com/watch?v=MsgpSl5JRbI

He is solving a LP by using a barrier function and has a very nice simulation
in the end where the solution converges to he optimal point and mu goes to zero.
I tried to replicate it in Maple but I am strugeling.
I have attached an worksheet what I have done so far:

InteriorPointMetho.mw

 
What I want to do is the same simulation in Maple as in the presentation for the below problem
ie how mu, objective function, and optimal allocation changes over time, starting from the global
optimal by I guess introducing slack variables. Any idea?!

restart:
with(plots):
ob := 100*x[1]*x[2]*exp((-x[1]^2)*(1/3)-(1/3)*x[2]^2):
con1 := 4*x[1]+5*x[2] <= 20:
con2 := 6*x[1]+x[2] <= 18:

x1 := evalf(solve(solve(convert(con1, equality), x[2]) = solve(convert(con2, equality), x[2]), x[1]), 5):
x2 := eval(solve(convert(con1, equality), x[2]), x[1] = x1):

my1 := plot(solve(convert(con1, equality), x[2]), x[1] = 0 .. 10, color = red, thickness = 2):
my2 := plot(solve(convert(con2, equality), x[2]), x[1] = 0 .. 10, color = green, thickness = 2):
my3 := contourplot(ob, x[1] = 0 .. 5, x[2] = 0 .. 5, color = black, contours = 10):
my4 := pointplot({[x1, x2]}, axes = boxed, symbol = solidcircle, symbolsize = 20):
my5 := textplot([x1, x2, "Optimal Allocation"], align = {ABOVE, RIGHT}, font = [times, bold, 16]):

[eval(ob, {x[1] = x1, x[2] = x2}), x[1] = x1, x[2] = x2];
display({my1, my2, my3, my4, my5}, view = [0 .. 4, 0 .. 5]);