(* ::Package:: *)

(*Coulomb.m*)

iota=10^(-8);(*This is a numerical trick used to prevent Vc from exploding at r = 0*)

Vc[i_,jprime,kprime]:=rhov*Deltax^3/(4*Pi*epsilon0*Sqrt[iota+(x[i]^2+z[kprime]^2+y[jprime]^2)]);
tVc=Table[{Sqrt[iota+(x[i]^2+z[kprime]^2+y[jprime]^2)],Vc[i,jprime,kprime]},{i,iini,ilast}];

(*Plot V(x) to check against analytical solution provided by Coulomb's law*)
(*t1D=Table[{Sqrt[x[i]^2+y[i]^2+z[i]^2],V[i,jprime,kprime]},{i,iini,ilast}];*)
t1D=Table[{Sqrt[iota+(x[i]^2+z[kprime]^2+y[jprime]^2)],V[i,jprime,kprime]},{i,iini,ilast}];

(*
Print[ListPlot[{t1D,tVc},PlotRange->{{-2,2},{0,1.5}},Joined->False,AspectRatio->1,AxesLabel->{"x(m)","V(x)"},
PlotStyle->{Black,Green},PlotLabel->"Numerical vs Analytical results (Black is for numerical result)"]];
*)
Print["For numerical result: Stepsize is ",Deltax, " Box size is  ", boxsize];
Print[ListPlot[{t1D,tVc},Joined->False,AspectRatio->1,AxesLabel->{"x(m)","V(x)"},
PlotStyle->{Black,Green},PlotLabel->"Numerical result in Black, Analytical in Green"]];
