Integrate[1/Sqrt[(x-x1)^2+(y-y1)^2],{x1,-1,1},{y1,-1,1},Assumptions->(x>1)&&(y>1)] 
f[x_, y_] = %
v = {D[-f[x, y], x], D[-f[x, y], y]}
ContourPlot[f[x, y], {x, .1, 2.1}, {y, .1, 2.1}, 
 ContourShading -> False, PlotPoints -> 50, AspectRatio -> Automatic,  Contours -> 20]
 VectorPlot[v, {x, .1, 2.1}, {y, .1, 2.1}]
 Show[%,%%]
 
 Plot[f[x,0],{x,.1,3}]
 Plot[4/x,{x,1.1,3}]
Show[%, %%, PlotRange -> All]

