My system consists of a coaxial conductor comprised of 2 concentric conductors. The outer conductor is earthed on its outer surface. The link shows the end of a semi-infinite coaxial such that the potential for most of it is given by that of an infinitely long coaxial conductor:
$\Phi(\phi,r,z) = \Phi_0 \frac{log(R_{outer}) - log(r)}{log(R_{outer})-log(R_{inner})} $
in cylindrical polar coordinates.
The Laplacian for this system seeing as though there is rotational symmetry becomes:
$\bigtriangledown^2 \phi = \frac{\partial^2 \phi}{\partial z^2}+\frac{1}{r} \frac{\partial \phi}{\partial r} +\frac{\partial^2 \phi}{\partial r^2} = 0$
and I have worked out that in difference notation this can be expressed in centred difference notation as follows (by using the same spatial resolution for z and r i.e. $\Delta z = \Delta r$):
$\phi_{i,k} =\frac{1}{4}(\phi_{i+1,k}+\phi_{i-1,k}+\phi_{i,k+1}+\phi_{i,k-1})+\frac{1}{8k}(\phi_{i,k+1}-\phi_{i,k-1})$ for r>0
and
$ \phi_{i,0}= \frac{2}{3} \phi_{i,1} + \frac{1}{6} \left(\phi_{i+1,0}+\phi_{i-1,0}\right)$ for r=0.
My question therefore is how do I arrange this into a form I can eventually solve on a computer with linear algebra methods, such as successive over-relaxation or Gauss-Siedel?

I misunderstood the Successive Over-Relaxation (SOR) method when asking this question.
To solve this an initial estimate of the value at a every point would have to be made.
The difference formulae can then be used to improve this estimate. Iterating over all points in space the formula can be used at each point using the initial values. These new values will converge on the correct answer given enough iterations. Once the value changes less than a arbitrarily defined tolerance from iteration to iteration it is accepted as the final answer (within an error given by the tolerance).