Thanks to all of you for giving a lot of knowledge and clear the confusion

Thank you Neil, I tried your suggestions above and worked well for up to the 3rd derivative. For the 4th derivative, the value of "ExactCondCovegeProbLAB22" is a negative large number which is not correct because it is a probability (it should be between 0 and 1). I plotted "funTable2D" and "ifun2D" to see how good the interpolation function "ifun2D" is as follows

ListPlot3D[funTable2Dc]

Plot3D[ifun2D[s, r], {s, LlimitxS, UlimitxS}, {r, LlimitxR, UlimitxR}]

As we can see, the two functions does not look similar

LlimitxR = height; UlimitxR = Sqrt[RmaxAB^2 + height^2];
LlimitxS = parmtrLAB[LlimitxR]; UlimitxS = parmtrLAB[UlimitxR];

Print[111];
funTable2D = 
  Table[{s, r, Exp[-noisepow*s]*laplaceLAB[s, r]}, {s, LlimitxS, 
    UlimitxS, 10000}, {r, LlimitxR, UlimitxR, 500}];
Print[222];
funTable2Dc = Flatten[Chop[funTable2D], 1];
Print[333];
ifun2D = Interpolation[funTable2Dc, Method -> "Spline"];
Print[444];
funDervtveLAB22[s_, r_, kk_] := (-s)^kk/kk!*D[ifun2D[s, r], {s, kk}];
xxSumS22[s_, r_] := 
  Sum[funDervtveLAB22[s, r, kk], {kk, 0, mLParm - 1}];
xxSumR22[r_] := xxSumS22[s, r] /. s -> (mLParm*sinr/zetaLAB*r^alfaLAB);
xxPdF22[r_] := xxSumR22[r]*CondPdfLAB[r];
ExactCondCovegeProbLAB22[[hIdx, mIdx]] = 
  NIntegrate[xxPdF22[r], {r, height, Sqrt[RmaxAB^2 + height^2]}];

Mohamed,

The problem you are having is that you are doing a numerical calculation and running out of precision. Numerical derivatives are noisy operations -- you must maintain extreme precision to take multiple derivatives. Your calculation involves a Numerical integrate and then you take the derivative of it. I recommend taking the derivative analytically -- it will also speed things up. Another way to do this is to keep precision with Interpolation functions:

fun = Table[{x, xxSumS1[x, 200]}, {x, 0, 500000, 1000}];

(This is a bit slow but you can lower the resolution. Also, You will get warnings which you can fix by specifying AccuracyGoal in the NIntegrate)

ifun = Interpolation[fun]
Plot[ifun[x], {x, 0, 500000}]
Plot[ifun'[x], {x, 0, 500000}]
Plot[ifun''[x], {x, 0, 500000}]

To get

Your calculations use NIntegrate, a numerical method. Symbolic differentiation cannot handle a numerical function.