Hi, I continue to have troubles while trying to use NDSolve for obtaining a very accurate solution (say, 20 accurate digits) of a large system of 240 first order ODEs. I have an independently obtained semi-analytical series solution for one unknown of this system, which should be valid for small values of the independent variable y. In fact, if I compare the numerical NDSolve solution with this series, I get a very good match (19-20 digits) for y > y0+delta and y< about 2, where y0=10^-30 is the initial y value, and delta (equal about 10^-28) is the first discrete integration step automatically selected by NDSolve. However, for y0 <y<y0+delta the accuracy is rather poor (no more than 10 digits). This looks as if there was some problem with the accuracy of the interpolation over the first integration step, within the Interpolating Function object returned by NDSolve. This happens despite requesting InterpolationOrder->All. My calling command is:

ivpsol=NDSolve[ivp,Table[Subscript[psi, i][y],{i,1,nnodes-1}],{y,y0,ymax},WorkingPrecision->50,PrecisionGoal->25,MaxSteps->Infinity,InterpolationOrder->All];

What could be the reason of such a behaviour and how to eliminate it? I tried to suppress the errors by adding StartingStepSize->10^-30 to the command, but this does not help, in fact this makes things worse, as the region of a poor accuracy then extends over several integration steps.

I need to integrate up to y=ymax=600, and all my other tests indicate that for such large y I obtain correct results. They appear accurate to at least 20 digits and converge to the asymptotic solution expected for large y. Hence, there is only a problematic issue with the first integration step.

Lesław