Ek = 0.04/1000;
PMASS = 0.938272;
gamm=(Ek+PMASS)/PMASS;
beta=Sqrt[1-(1/gamm^2)];
dT = 10 *10^-12;
q:=1.6021773349*10^-19
m:=1.6726231*10^-27;
Q:=13.4*10^-12;
c:=2.99792 10^8;
E0:=8.85418782*10^-12;
maxpos:=1;
maxtime:=maxpos/(beta*c)
k:=(3*q*Q)/(4*Pi*m*gamm^3*E0);
myScale:=1.0;

Mx[x_,y_,z_]:=0.5*NIntegrate[x*y*gamm*z/((x^2+s)*((x^2+s)*(y^2+s)*(gamm^2*z^2+s))^0.5),{s,0,Infinity}];
My[x_,y_,z_]:=0.5*NIntegrate[x*y*gamm*z/((y^2+s)*((x^2+s)*(y^2+s)*(gamm^2*z^2+s))^0.5),{s,0,Infinity}];
Mz[x_,y_,z_]:=0.5*NIntegrate[x*y*gamm*z/((gamm^2*z^2+s)*((x^2+s)*(y^2+s)*(gamm^2*z^2+s))^0.5),{s,0,Infinity}];

solution=NDSolve[{x''[t]-k*Mx[x[t],y[t],z[t]]/(y[t]*z[t])==0,x[0]==0.001,x'[0]==0,y''[t]-k*My[x[t],y[t],z[t]]/(x[t]*z[t])==0,y[0]==0.001,y'[0]==0,z''[t]-k*Mz[x[t],y[t],z[t]]/(x[t]*y[t])==0,z[0]==0.0001,z'[0]==0},{x,y,z},{t,0,maxtime},Method->Automatic];

LogPlot[{
Evaluate[(1/(Sqrt[5]))*x[t*maxtime/maxpos]/.solution],
Evaluate[(1/(Sqrt[5]))*y[t*maxtime/maxpos]/.solution],
Evaluate[(1/(Sqrt[5]))*z[t*maxtime/maxpos]/.solution]},
{t,0,maxpos},
PlotPoints->300,
PlotRange->{0,0.03},
PlotLabel->"Rms Analytical Solution d=1m, Q=13.4 pC",
GridLines->Automatic,
AxesLabel->{"z(m)","R(m)"}
]
(*Export["round870keVRms.eps",%]*)


mydata1=Table[ Flatten[
N[{ t, Evaluate[(1/(Sqrt[5]))*x[myScale*maxtime*t/maxpos]/.solution],
          Evaluate[(1/(Sqrt[5]))*y[myScale*maxtime*t/maxpos]/.solution],
          Evaluate[(1/(Sqrt[5]))*z[myScale*maxtime*t/maxpos]/.solution]}]],{t,0,maxpos,maxpos/3000}];
(* Export["round870keVRms.dat",mydata1,"Table"]*)