The code used on this page is availible here.
I was kindly sent some unpublished code (related to Möbius function) by Franz Lehner with permission to merge this with my own poset related code described here.
These are my notes related to the code written by Franz.
The Möbius function can be applied to any type that factorises
As an example take the factors of 12. An arrow means 'has factor'. I have drawn two arrows from 4 to 2 just to emphasise that 4 is 2 squared, but of course a poset has, at most, one arrow between nodes. 
So here is the code to generate a poset for this:
(1) > FP := FinitePoset(NNI) (1) FinitePoset(NonNegativeInteger) Type: Type (2) > carrier:List(NNI) := [1,2,3,4,6,12] (2) [1,2,3,4,6,12] Type: List(NonNegativeInteger) (3) > (x:NNI,y:NNI):Boolean +> ((x rem y) < 1) Function declaration isaFactor : (NonNegativeInteger, NonNegativeInteger) > Boolean has been added to workspace. Type: Void (4) > pset1 := finitePoset(carrier,isaFactor)$FP Compiling function isaFactor with type (NonNegativeInteger, NonNegativeInteger) > Boolean +true true true true true true+   false true false true true true   false false true false true true (4)   false false false true false true   false false false false true true   +false false false false false true+ Type: FinitePoset(NonNegativeInteger) 
Incidence Algebra
At the moment my (Martin Baker) poset code uses an adjacency matrix for the Rep of a poset (with entries being boolean values  so an adjacency array). See FinitePoset code here.
For the Incidence Algebra the entries can be any ring such as integer. For Zeta function the values are limited to 0 and 1 which can be easily mapped to false and true.
Franz Lehner originally used the name IncidenceMatrix for this. I changed the name because it tends to be used for matrix where: rows are vertices and columns are edges as described here. So I changed the name IncidenceMatrix to IncidenceAlgebra.
An incidence algebra is an algebra over a field with convolution as its product. This can be defined in terms of a finite poset as follows:
If P is a finite poset then the incidence algebra is:
I(P) = { α : P × P> α(x,y) = 0 if not x ≤ y}
This has the following operations:
 Addition: pointwise addition: (α + β)(x,y) = α(x,y) + β(x,y)
 Scalar multipication: (kα)(x,y) = k•α(x,y) for k
 Convolution: (α * β)(x,y) = Σ_{z}α(x,z) β(z,y)
So addition and scalar multipication are the same as the axioms for vector algebra and convolution is similar to matrix product.
Zeta Function
An entry of 1 means relation is true and 0 means it is not.
ζ(x,y) = {  0 if not x ≤ y  } 
1 if x ≤ y 
The zeta function produces a matrix [multipication table with only 0 and 1]
(5) > z := zetaMatrix(pset1) +1 1 1 1 1 1+   0 1 0 1 1 1   0 0 1 0 1 1 (5) [ ,[1,2,3,4,6,12]] 0 0 0 1 0 1   0 0 0 0 1 1   +0 0 0 0 0 1+ Type: IncidenceAlgebra(Integer,NonNegativeInteger) 
Now try generating zeta matrix again but this time we put carrier set in a different order:
(6) > jumble:List(NNI) := [12,2,3,1,6,4] (6) [12,2,3,1,6,4] Type: List(NonNegativeInteger) (7) > pset2 := finitePoset(jumble,isaFactor)$FP +true false false false false false+   true true false false true true    true false true false true false (7)   true true true true true true    true false false false true false   +true false false false false true + Type: FinitePoset(NonNegativeInteger) (8) > zetaMatrix(pset2) +1 1 1 1 1 1+   0 1 0 1 0 1   0 0 1 1 1 1 (8) [ ,[1,3,2,6,4,12]] 0 0 0 1 0 1   0 0 0 0 1 1   +0 0 0 0 0 1+ Type: IncidenceAlgebra(Integer,NonNegativeInteger) 
We can see that our carrier set has been reordered. It is not exactly in the same order as before, 2,3 have been reversed and 4,6 have been reversed. However both of these are compatable with the topological ordering because we could take either the left or the right branch on the diagram first. 
Möbius Function
As explained on Wikipedia here, the value of the Möbius function is −1, 0 or 1 depending on the factorization of n into prime factors:
The first few values are given in the table here: So, how does this relate to the matrix? 

μ = ζ^{1}
The Möbius function is the inverse of the zeta function.
This gives the axiom:
δ(x,y) = (μ * ζ)(x,y) = Σ_{z∈(x,y)}μ(x,z) ζ(z,y)= Σ_{z∈(x,y)}μ(x,z)
μ(y) = {  1  if y = 0  } 
Σ_{z<y}μ(z)  if y < 0 
Now generate Möbius matrix:
(9) > m := moebius(pset1) +1  1  1 0 1 0 +   0 1 0  1  1 1    0 0 1 0  1 0  (9) [ ,[1,2,3,4,6,12]] 0 0 0 1 0  1   0 0 0 0 1  1   +0 0 0 0 0 1 + Type: IncidenceAlgebra(Integer,NonNegativeInteger) (10) > z*m +1 0 0 0 0 0+   0 1 0 0 0 0   0 0 1 0 0 0 (10) [ ,[1,2,3,4,6,12]] 0 0 0 1 0 0   0 0 0 0 1 0   +0 0 0 0 0 1+ Type: IncidenceAlgebra(Integer,NonNegativeInteger) 
As we can see in (line 10) the zeta multipied by the Möbius matrix produces the unit (delta) matrix.
Cover Matrix
Returns all elements covered by x arithmetics
The covering relation is determined as follows.
 We assume that the list of elements xx is topologically sorted.
 We do a double loop over x in xx and y in xx.
Since the possible covers of an element come later in the list, it suffices to let y run over the rest of the list. Moreover, any y ≥ x is a cover unless some predecessor lies in between.
Now generate cover matrix:
(11) > coverMatrix(pset1) +0 1 1 0 0 0+   0 0 0 1 1 0   0 0 0 0 1 0 (11) [ ,[1,2,3,4,6,12]] 0 0 0 0 0 1   0 0 0 0 0 1   +0 0 0 0 0 0+ Type: IncidenceAlgebra(Integer,NonNegativeInteger) 