(* ::Package:: *)

(*simpson.m*)
(*Eq. E.11*)

(*Define Odd NOdd*)
NOdd=101;z[0]=-L;z[NOdd]=L;

(*define the integral*)
tBysimp={};

For[i=0,i<=NStepx,i++,
x = xini + i*Deltax;
sumfodd=0;sumfeven=0;


For[k=0,k<=(NOdd-1)/2,k++,
indexodd=2*k+1;
indexeven=2*(k-2)+4;
zodd=z[0]+indexodd*Deltaz;
zeven=z[0]+indexeven*Deltaz;
fodd=(mu0I/(4*Pi))*x/(zodd^2+x^2)^(3/2);
feven=(mu0I/(4*Pi))*x/(zeven^2+x^2)^(3/2);
sumfodd=sumfodd+fodd;
sumfeven=sumfeven+feven;
(*Print["indexodd ",indexodd," indexeven, ",indexeven, " k ",k];*)
];(*end for k*)


sumfodd=(4*Deltaz/3.0)*sumfodd;
sumfeven=(2*Deltaz/3.0)*sumfeven;
sumendpoints=(2*Deltaz/3.0)*(f[x,z[0]]+f[x,z[NOdd]]);
sumBy=sumfodd+sumfeven+sumendpoints;

tBysimp=Append[tBysimp,{x,sumBy}];
(*Print["x=",x];
Print[tBysimp];*)
](*end for i*)
