The following is a simple code (for solving just one ODE), which illustrates my approach:

wprec=50;
tmax1=1;
tmax2=10;
ivpsol1=NDSolve[{u'[t]+u[t]^2==0,u[0]==1},{u[t]},{t,0,tmax1},WorkingPrecision->wprec,PrecisionGoal->wprec/2,
MaxSteps->Infinity,InterpolationOrder->All];
unum1[t_]:=Evaluate[First[u[t]/.ivpsol1]];
umax1=unum1[tmax1]
0.49999999999999999999999936808923840799082795395903
Export["Z:\\umax1.txt",umax1,"List"];
umax1new=First[SetPrecision[Import["Z:\\umax1.txt","List"],wprec]]
0.49999999999999999999999936808923840799082795395903
ivpsol2=NDSolve[{u'[t]+u[t]^2==0,u[tmax1]==umax1new},{u[t]},{t,tmax1,tmax2},WorkingPrecision->wprec,PrecisionGoal->wprec/2,
MaxSteps->Infinity,InterpolationOrder->All];
unum2[t_]:=Evaluate[First[u[t]/.ivpsol2]];
unum2[tmax1]
0.4999999999999999999999993680892384080

In this simple case I don't get warnings indicating that the precision of my ODE is smaller than the WorkingPrecision during the second call of NDSolve. It seems also that there is no problem with preserving all 50 significant digits of umax1new during Export/Import. However, the Interpolating Function object resulting from the second NDSolve call does not reproduce all digits of the initial condition (umax1new), although more than 25 requested digits are reproduced.

By the way, I observe that for this example NDSolve successfully maintains the requested precision of the solution for t up to about 1 (exact value of umax1new is 0.5) , but later the solution obtained becomes less and less accurate, independently of whether I restart the calculations or not. The adaptive step selection mechanism does not seem to work properly for larger t.

Leslaw

From your description, I assume it's loss of precision from evaluating the boundary conditions to be used as initial conditions for the next stage. This loss comes from an estimate of the maximum roundoff error propagated through the computation. It is an estimated bound, not the actual error.

Setting $MinPrecision and $MaxPrecision should suppress the message. If your BCs are numerically ill-conditioned, it won't help that; it will probably hide it (caveat).

Block[{$MinPrecision = 50, $MaxPrecision = 50}, BCs /. currentstate]

(There seem to be two, distinct kinds of precision being discussed, namely, the floating-point precision of the numbers used by the kernel and the precision (or accuracy) of the solution. Whenever I do what is described, the solution at the end of a stage always has a floating-point precision of 50, no matter what the precision of the solution happens to be. In fact, NDSolve[] cannot know the precision of the solution. It knows only the format of the numbers (50-digit bignums) in its argument. So NDSolve[] is not complaining about the solution being inaccurate. It's complaining about some of the numbers having a bignum format of less than 50 digits.)