phi[x_,z_]=-z(1+1/Sqrt[x^2+z^2]^3)

v={-D[phi[x,z],x],-D[phi[x,z],z]}

StreamPlot[v,{x,-2,2},{z,-2,2},RegionFunction->Function[{x,z},x^2+z^2>1.3]]

RowReduce[{{5,3,18},{3,9,6}}]

RowReduce[{{3,0,2,14},{-1,2,1,0},{0,3,-1,14}}]