If we given two parametric lines, defined as l & k (with points
lStart, lEnd, kStart, kEnd); how do we find the shortest line (j,
defined by jStart, jEnd) joining two infinite skew lines?

I'd approach it like this:

Generic point on line l: x(t) = lStart + t (lEnd - lStart)
Generic point on line k: y(s) = kStart + s (kEnd - kStart)

Then let f(s, t) = ||(x(t)-y(s))||^2 = squared length of line j
and minimise this w.r.t. s and t.

E.g. let lStart = (1, 0, 1), lEnd = (2, 1, 1), kStart = (-1, 1, -2)
and kEnd = (0, 0, 1). Then we get
x(t) = (1, 0, 1) + t (1, 1, 0) = (1 + t, t, 1)
y(s) = (-1, 1, -2) + s (1, -1, 3) = (s - 1, 1 - s, 3s - 2)

f(s, t) = (2+t-s)^2 + (t+s-1)^2 + (3s-3)^2

(Note: Here I will use df/ds, df/dt to refer to partial derivatives,
not total ones as this notation normally signifies - due to the lack
of a delta symbol in ASCII!)

df/dt = 2(2+t-s).1 + 2(t+s-1).1
df/ds = 2(2+t-s).(-1) + 2(t+s-1).1 + 2(3s-3).3

For a minimum (and it always will be a minimum) we need both of these
to be 0. Simplifying, we get:

{  4t + 2  = 0
{ 22s - 24 = 0

giving t = -1/2, s = 12/11. Therefore the points jStart and jEnd are given by

x(t) = (1 + t, t, 1) = (1/2, -1/2, 1)
y(s) = (s - 1, 1 - s, 3s - 2) = (1/11, -1/11, 14/11).

You could probably go through analytically and get an explicit formula
for this if you really wanted to - let me know if you want me to have
a bash at it. Though someone's almost certainly put a formula up on
the web somewhere.

I agree with your solution, but I was hoping for more of a geometric
approach using a solution that could be established using projections
/ converse projects (without derivatives). I'm quite aware of how to
use projects / converse projections, however the material which I have
available to me is lacking.