equation solving – NDSolve with embedded FindRoot and NIntegrate


I am trying to plot results for a very simple differential equation of the form:

$$frac{partial^2 x(N,z'(N))}{partial N^2} = F(N,z'(N)), $$

where $z'(N)$ is a function of $N$ that needs to be solved using FindRoot for at every $N$ position, and $F(N,z’)$ is a nasty equation that results from a numerical integration over:

$$ F(N,z’) = int_{-infty}^{infty} expleft( -frac{x’^2}{2sigma_{x’}^2} right) F(N,z’,x’)dx’$$.

So, I put together some mathematica code but it runs terrible slow (on the order of a day or two)! I noticed there were some things that affected the speed of the code particularly the numerical coefficient in front of $F(N,z'(N))$. But I was wondering if there is any help to be given to get better/faster results! Any help would be greatly appreciated!

Note: I had to use $NN$ in place of $N$ in my equations because its a function in mathematica. Also, in the FN function I have to actually copy and paste the output of FNzprime (an ugly mess) in the integrand for it to evaluated.

(*constants*)
e = -1.60217733*10^-19;
m = 9.109389699999999*10^-31;
epsilon = 8.854187817620391*10^-12;

(*basic equations*)
rs2 = {zprime, xprime + K/(gamma*kw) Sin(kw*zprime), 0};
ro2 = {(NN + 10000)*lw, x + K/(gamma*kw) Sin(kw*(NN + 10000)*lw), 0};

betas = {beta - K^2/(4 gamma^2) Cos(2 kw*zprime),K/gamma Cos(kw*zprime), 0};
betao = {beta - K^2/(4 gamma^2) Cos(2 kw*(NN + 10000)*lw),K/gamma Cos(kw*(NN + 10000)*lw), 0};

betaDot = {(c*K^2*kw)/(2 gamma^2)Sin(2 kw*zprime), -((c*K*kw)/gamma) Sin(kw*zprime), 0};

deltar2 = ro2 - rs2;
Rgam2 = Sqrt(deltar2((1))^2 + deltar2((2))^2);
Ec2 = (e/(4 Pi*epsilon)) (deltar2/Rgam2 - betas)/(gamma^2 Rgam2^2 (1 - (deltar2/Rgam2).betas)^3);
Erad2 = (e/(4 Pi*epsilon)) Cross(deltar2/Rgam2, Cross(deltar2/Rgam2 - betas, betaDot))/(c*Rgam2*(1 - (deltar2/Rgam2).betas)^3);

Bc2 = Cross(deltar2/Rgam2, Ec2);
Brad2 = Cross(deltar2/Rgam2, Erad2);

Fbc2 = Cross(betao, Bc2);
Fbrad2 = Cross(betao, Brad2);


sumEtran = (Ec2((2)) + Erad2((2)));
sumFBtran = Fbc2((2)) + Fbrad2((2));



(*Numeric Functions*)

ZPRIME(NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, gamma_, K_, kw_, beta_, sigma_, lw_) :=zprime /. FindRoot(sigma == (1/(gamma kw))Sqrt(gamma^2 + K^2) (EllipticE(kw*(NN + 10000)*lw, K^2/(gamma^2 + K^2)) - EllipticE(kw zprime, K^2/(gamma^2 + K^2))) - beta (Sqrt)(((NN + 10000)*lw - zprime)^2 + (x - xprime + (K Sin(kw *(NN + 10000)*lw))/(gamma kw) - (K Sin(kw zprime))/(gamma kw))^2), {zprime, 0})


coeff = ((e*lw^2)/(gamma*m*beta^2*c^2) (10^-10/e)/(2 Pi (30*10^-6) (10^-5)) Exp(-(sigma^2/(2 (10^-5)^2))));
FNzprime =coeff (sumEtran + sumFBtran) /. {lw -> 0.026, K -> 1, beta -> Sqrt(1 - 1/(4000/0.511)^2), gamma -> 4000/0.511, c -> 3*10^8, kw -> (2 Pi)/0.026, zprime -> ZPRIME}

FN(NN_?NumericQ, x_?NumericQ, sigma_?NumericQ) :=With({ZPRIME = ZPRIME(NN, x, 0, 4000/0.511, 1, (2 Pi)/0.026, Sqrt(1 - 1/(4000/0.511)^2), sigma, 0.026)}, 
  NIntegrate( (Exp(-(xprime^2/(2 (30*10^-6)^2)))) FNzprime, {xprime, -300*10^-6, 300*10^-6}))

sol00 = NDSolve({X''(NN) - (FN(NN, 0, 10^-8)) == 0, X(0) == 0, X'(0) == 0}, X, {NN, 0, 140})

Plot(X(NN) /. {sol00}, {NN, 0, 10}, Evaluated -> True)