
phi=-1/Sqrt[1+z^2]
dphi=D[phi,z]
ddphi=D[dphi,z]
            
Plot[phi,{z,-10,10},AxesLabel->{"z","phi"},PlotRange->{{-10,10},{0,-1}}]
Plot[-dphi,{z,-10,10},AxesLabel->{"z","Ez"}]
Plot[ddphi,{z,-10,10},AxesLabel->{"z","Er"},PlotRange->Full]

z0=100
r0=0.01
equation={4*phi*r''[z]+2*dphi*r'[z]+ddphi*r[z]==0,r[-z0]==r0,r'[-z0]==0}
sol = NDSolve[equation,{r[z]},{z,-z0,z0}]
            
Plot[r[z]/.sol,{z,-z0,z0},AxesLabel->{"z/a","r/a"}]
Plot[r[z]/.sol,{z,-10,10},AxesLabel->{"z/a","r/A"}]
