How to pick out specific parts of the solution of a differential equation efficiently?

Question

I'm trying to solve a set of differential equations (in parallel) and then select specific parts of the solution. What I have at the moment (see below-I've simplified the differential equations, in fact they include matrices so that xx, yy and zz will be matrices) works, but uses a lot of memory. Ultimately, I only need the solutions of the differential equations at particular points- is there an efficient way to do this so that Mathematica only remembers the parts I need?

Thanks!

X = SparseArray[KroneckerProduct[{{0, 0, 0, 0, 0}, {0.5, 0, 0, 0, 0}, {1, 0, 0, 0,0}, {0, 0, 0, 0, 0}, {0, 1, 0.5, 0, 0}}, IdentityMatrix[3]]];
\[Rho]g = KroneckerProduct[{{0 , 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 1}}, {{1, 0, 0}, {0, 0, 0}, {0,
  0, 0}}];
sol = ParallelTable[NDSolve[{xx'[t] == xx[t], yy'[t] == c*yy[t], zz'[t] == zz[t],xx[0]==\[Rho]g, yy[0] == \[Rho]g, zz[0] == \[Rho]g}, {xx, yy, zz}, {t, 0, 2000}], {c, 0, 1000,10}];
a1[x_, i_] := (xx[t] /. sol[[i, 1]])
a2[x_, i_] := (yy[t] /. sol[[i, 1]])
a3[x_, i_] := (zz[t] /. sol[[i, 1]])
answer[t1_, i_] := Total[X.Transpose[a1[t1, i] - a2[t1, i] - a3[t1, i]], 2]
choosevalues = Table[answer[t1, i], {i, 1, 101}, {t1, 500 + (i - 1)*10, 500 + (i -1)*10 + 1000, 10}];

Answer 1

Some food for brain! Will elaborate soon.

Clear[sol];
sol[tval_?NumericQ, cval_?NumericQ] := 
Module[{res, c}, 
  res = ParametricNDSolveValue[{xx'[t] == xx[t], yy'[t] == c*yy[t], 
  zz'[t] == zz[t], xx[0] == 2, yy[0] == 1, 
  zz[0] == 2}, {xx[tval] - yy[tval] - zz[tval]}, {t, 0, 100}, {c},
  MaxSteps -> 10^6, AccuracyGoal -> 2];
  res[cval]
  ];

This plays the role of your answer[t,i]

sol[15, 12]

{-1.53353*10^78}

Update:

First lets handle the equations properly. I inserted the HilbertMatrix to get nontrivial results. Your example matrix is too full of zeros.

var={xx[t],yy[t],zz[t]}=Table[Array[i[#1,#2][t]&,{15,15}],{i,{x,y,z}}];
{xx'[t],yy'[t],zz'[t]}=({xx[t],yy[t],zz[t]})/.a___[t]->a'[t]; 
X=SparseArray[KroneckerProduct[{{0,0,0,0,0},{0.5,0,0,0,0},{1,0,0,0,0},{0,0,0,0,0},  
 {0,1,0.5,0,0}},IdentityMatrix[3]]];
ρg=KroneckerProduct[HilbertMatrix[5],{{1,0,0},{0,0,0},{0,0,0}}];
A = RandomReal[{-1, 1}, {15, 15}];
B = c RandomReal[{-.5, .5}, {15, 15}]; (* B Depends on c *)
Com[x_, y_] := x.y - y.x;
eqs=Flatten@(Join@@MapThread[#1==#2&,#,2]&/@{
   (* Equations {LHS matrix, RHS matrix} *)
   (* the xx'[t]==Com[(A-B[c]),xx[t]] *)
   {xx'[t],Com[(A - B), xx[t]]},
   {yy'[t],c yy[t]},
   {zz'[t],zz[t]},
   (* Equations LHS\[Equal]RHS with {LHS matrix & RHS matrix} *)
   {xx[t],ρg}/.t-> 0,
   {yy[t],ρg}/.t-> 0,
   {zz[t],ρg}/.t-> 0});

Now the version 8 compatible function sol8 that takes different values of c. The process can be made more optimal. Look at the suggested post by @Kuba.

Clear[sol8];
sol8[tval_?NumericQ, cval_?NumericQ] := Module[{rule},
   rule =Dispatch@Flatten@NDSolve[
         Evaluate[eqs /. c -> cval], 
         Flatten@var /. a___[t] -> a,
         {t, 0, 20},AccuracyGoal -> 2];
  (X . (xx[t] - yy[t] - zz[t])) /. rule /.t -> tval]; 

Now testing!

DistributeDefinitions[sol8];
choosevalues = ParallelTable[Evaluate@sol8[tval, -Cos@cval], 
 {cval, 1, 10, 1}, {tval, 1, 20, .5}];

Check the dimension of choosevalues

Dimensions@choosevalues 

{10, 39, 15, 15}

In Mathematica 9 you can use Image3D to visualize the solution matrices as voxels.

Image3D[choosevalues]

Hope this gets you going!

BR

Answer 2

Yes, there is an efficient way. It is important to see that the range of solution does not have to start from where initial conditions are. So you can request a solution to be between any points, far removed from the initial conditions. Here is an example:

ic = y[0] == 1;
ode = y'[x] == y[x] Cos[x + y[x]];
sol = NDSolve[{ode, ic}, y, {x, 0, 30}]
Plot[y[x] /. sol, {x, 0, 30}, PlotRange -> All, Evaluated -> True]

But suppose you want the solution only from x=15 to x=16 ? You can just ask for that part, and still use the same initial conditions at x=0

ic = y[0] == 1;
ode = y'[x] == y[x] Cos[x + y[x]];
sol = NDSolve[{ode, ic}, y, {x, 15, 16}]
Plot[Evaluate[y[x] /. sol], {x, 15, 16},PlotRange -> {{0, 30}, {0, 1}}, PlotStyle -> Red]

You can even do more. You can as NDSolve to solve for small segments, and treat it as an integration step itself, and use the new solution to solve for the next segment.