I am trying to use NDSolve to solve some a set ODEs coupled in a slightly slightly unusual, and am running into some problems.

The following code generates the equations of interest, where one the solution of one ODE is used as the initial condition for a second, generating a collection of 2D interpolating functions (apologies for the large matrix, this is the minimal example I could think of!):

(*Prepare vector of dependent variables*)
Rho[t_]= Flatten@Array[Subscript[r,#1,#2][t]&,{3,3}];
Lamb[t_, tau_]=Flatten@Array[Subscript[La,#1,#2][t,tau]&,{3,3}];

(*Include propagator matrix*)
lio ={{-0.0000506308, 0. + 1. I, 0, 0. - 1. I, 0, 0, 0, 0, 0}, {0. + 1.I, -9.72405, 0, 0, 0. - 1. I, 0, 0, 0, 0}, {0,  0, -4.72405, 0, 0, 0. - 1. I, 0, 0, 0}, {0. - 1. I, 0, 0, -9.72405,   0. + 1. I, 0, 0, 0, 0}, {0, 0. - 1. I, 0, 0. + 1. I, -10., 0, 0, 0,   0}, {0, 0, 0. - 1. I, 0, 0, -5., 0, 0, 0}, {0, 0, 0, 0, 0,   0, -4.72405, 0. + 1. I, 0}, {0, 0, 0, 0, 0, 0, 0. + 1. I, -5.,   0}, {0.0000506308, 0, 0, 0, 10., 0, 0, 0, 0}};

(*Solve the first ODE*) 
rho0 = {1, 0, 0, 0, 0, 0, 0, 0, 0};
Rsol = First@NDSolve[{D[Rho[t],t]==lio.Rho[t], Rho[0]==rho0}, Rho[t],{t,0,100}];

(*Solve the second ODE using the first as initial condition*) 
Lam0 = {0., 0., 0., 0., 0., 0., 0. + 1. Subscript[r, 2, 1][t], 0. + 1. Subscript[r, 2, 2][t], 0. + 1. Subscript[r, 2, 3][t]}/.Rsol;
LambSol = NDSolveValue[{D[Lamb[t,tau],tau]==lio.Lamb[t,tau], Lamb[t,0] == Lam0},Lamb[t,tau],{t,0,100},{tau,0,3}];

The issue arises when we compare the solutions in the t variable. Plotting RSol[[5]] and the matching element LambSol[[-2]], with tau = 0, we find:

Show[LogLinearPlot[{Evaluate[Rho[t][[5]]/.Rsol],Evaluate[LambSol[[-2]]/.tau->0]},{t,10^-3,20}],
ListLogLinearPlot[Transpose[{LambSol[[-2,0]]["Coordinates"][[1]],
LambSol[[-2,0]]["ValuesOnGrid"][[All,1]]}],PlotStyle->Red]]

As you can see the dynamics found from the second ODE solution (yellow) doesn’t match the converged solution from RSol (blue). I think the reason is that NDSolve is not taking enough points the along the t-axis in second NDSolve (given by the red points). I can fix this to an extent by decreasing MaxStepFraction, however taking this too small or for systems of equations that are too large makes the code take excessive amounts of time, or the Kernel crash.

Does anybody have insight into using NDSolve in this way? I’m at a loss of how to tackle this problem.