sol=Solve[{a1==d1+e1,d1+e e1==c1-h0,a1==k (d1-2 e1),k (d1-2 e e1)==-2 c1-h0},{a1,d1,e1,c1}]

sol1=First[sol]/.{e->.8^3,h0->1,k->100}
sol1=Join[sol1,{e->.8^3,h0->1,k->100}]

v[x_,z_]:=If[x^2+z^2<e^(2/3),a1 z,If[x^2+z^2<1,d1 z+e e1 z/Sqrt[x^2+z^2]^3,c1 z/Sqrt[x^2+z^2]^3-h0 z]]/.sol1

Plot[v[0,z],{z,-2,2}]
ContourPlot[v[x,z],{x,-2,2},{z,-2,2},Contours->15,ContourShading->False]

b1[x_,z_]=Evaluate[{D[d1 z+e e1 z/Sqrt[x^2+z^2]^3,x],D[d1 z+e e1 z/Sqrt[x^2+z^2]^3,z]}]/.sol1
b2[x_,z_]=Evaluate[-{D[c1 z/Sqrt[x^2+z^2]^3,x],D[c1 z/Sqrt[x^2+z^2]^3,z]-h0}]/.sol1

b[x_,z_]:=If[x^2+z^2<e^(2/3),{0,-a1},If[x^2+z^2<1,-k b1[x,z],b2[x,z]]]/.sol1
VectorPlot[b[x,z],{x,-2,2},{z,-2,2}]