Here is my code:

k = 0.42; e = 5; d = 7; n = 1/0.7158; R2 = 40000;

s =
NDSolve[{k^2  (u'[y] )^(3/2) (u'[y] + 2 u[y] Exp[y]/(R2 - Exp[y]))^(
1/2) - (-u'[y] (R2 Exp[-y] - 1) - 2 u[y]) (n - 1)/R2 - 1 == 0,
u[Log[4.17 (n - 1)]] == 4.17},
u, {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}];

C1[y_] = Evaluate[u[y] /. s[[2]]];

Plot[C1[y], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

x = (((-C1'[y] (R2 Exp[-y] - 1) - 2 C1[y])/(-Exp[y] + R2)^3 +
R2/((n - 1) (-Exp[y] + R2)^3)) Abs[
C1'[y] Exp[-y] + 2 C1[y]/(R2 - Exp[y])])^(-1/4) Exp[-y] (n - 1)^(
1/2) R2^(-1/2)

f[y_] = NIntegrate[
2/3 a^(-5/3)
Exp[-a d Abs[x]] (1 + (e n Exp[y]/(R2 (n - 1)))^2/a^2)^(-17/
6), {a, 1, Infinity}]

f1 = FunctionInterpolation[(f[y])^(
3/4), {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

Plot[Re[f1[y]], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

s = NDSolve[{k^2  (f1[y]) (u'[y])^(
3/2) (u'[y] + 2 u[y] Exp[y]/(R2 - Exp[y]))^(
1/2) - (-u'[y] (R2 Exp[-y] - 1) - 2 u[y]) (n - 1)/R2 - 1 == 0,
u[Log[4.17 (n - 1)]] == 4.17},
u, {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

C1[y_] = Evaluate[u[y] /. s[[2]]]

Plot[Re[C1[y]], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

x = (((-C1'[y] (R2 Exp[-y] - 1) - 2 C1[y])/(-Exp[y] + R2)^3 +
R2/((n - 1) (-Exp[y] + R2)^3)) Abs[
C1'[y] Exp[-y] + 2 C1[y]/(R2 - Exp[y])])^(-1/4) Exp[-y] (n - 1)^(
1/2) R2^(-1/2)

f[y_] = NIntegrate[
2/3 a^(-5/3)
Exp[-a d Abs[x]] (1 + (e n Exp[y]/(R2 (n - 1)))^2/a^2)^(-17/
6), {a, 1, Infinity}]

f1 = FunctionInterpolation[(f[y])^(
3/4), {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

Plot[Re[f1[y]], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

(*up to here, I think there is no problem. then the kernel crashes below*)

s = NDSolve[{k^2  (f1[y]) (u'[y])^(
3/2) (u'[y] + 2 u[y] Exp[y]/(R2 - Exp[y]))^(
1/2) - (-u'[y] (R2 Exp[-y] - 1) - 2 u[y]) (n - 1)/R2 - 1 == 0,
u[Log[4.17 (n - 1)]] == 4.17},
u, {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

C1[y_] = Evaluate[u[y] /. s[[2]]]

Plot[Re[C1[y]], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]


R2 is a constant. When I set R2 smaller than 33000, I can always got a numerical solution of u[y]. However when R2 gets very large,i.e. 40000, mathematica just doesn't work at all. I mean there is no solution given for the third differential equation.

Can anybody tell me how to solve this problem? I will be more than appreciated.

Thanks.

-
Please post a fully self contained non-working example that one can copy/paste and run to show the problem including any numerical values used. –  Nasser Mar 11 '14 at 0:57
Thanks for your reply. I have posted the full code. –  zeo Mar 11 '14 at 1:07
Can you reduce the size of your example, and indicate where exactly and what goes wrong? It'll save people a lot of time and might get you an answer. All I see now is that at some point the kernel crashes, but I didn't take the time to find out at which step. I think you can do this yourself and let us know. –  Szabolcs Mar 11 '14 at 1:32
That said, a kernel crash indicates a bug, so can you report this to support at wolfram.com ? –  Szabolcs Mar 11 '14 at 1:34
Thanks, I edited the post again, indicating where exactly goes wrong –  zeo Mar 11 '14 at 1:45