Как решить систему интегральных уравнений типа "свертка"?

В случае аппроксимации свертки методом прямоугольников получаю производящую функцию системы линейных уравнений (в блокноте представлена матрицей А).

Шаг суммирования 'Nu' с заданным интервалом h. Получаю матрицу, размером 7 х n, где 7 есть количество уравнений в системе, а Nu,-количество столбцов (количество систем, эволюция которых зависит от шага). Аналитически решил в образах Лапласа в разных вариантах, а вот численно не получается, так как не могу зациклить решатель Solve в смысле подстановки значений начального решения (1-й столбец) во второй и так далее рекурсивно. Промежуточные матрицы p и fun формируют ядра уравнений и заданы. Начальные условия {1,0,0,0,0,0,0}.

h := 0.01;
f1[\[Lambda]1_, t_] := \[Lambda]1*Exp[-\[Lambda]1*t];
f2[\[Lambda]2_, t_] := \[Lambda]2*Exp[-\[Lambda]2*t];
f3[\[Lambda]3_, t_] := \[Lambda]3*Exp[-\[Lambda]3*t];
f4[\[Lambda]4_, t_] := \[Lambda]4*Exp[-\[Lambda]4*t];
f5[\[Lambda]5_, t_] := \[Lambda]5*Exp[-\[Lambda]5*t];
f6[\[Lambda]6_, t_] := \[Lambda]6*Exp[-\[Lambda]6*t];
f7[\[Lambda]7_, t_] := \[Lambda]7*Exp[-\[Lambda]7*t];
f8[\[Lambda]8_, t_] := \[Lambda]8*Exp[-\[Lambda]8*t];
f9[\[Lambda]9_, t_] := \[Lambda]9*Exp[-\[Lambda]9*t];
f10[\[Lambda]10_, t_] := \[Lambda]10*Exp[-\[Lambda]10*t];
p = ( {
    {1, 1, 0, 0, 0, 0, 0},
    {0, 0, 1, 0, 0, 0, 0},
    {0, 0.3, 0, 0.6, 0.1, 0, 0},
    {0, 0.4, 0, 0, 0.1, 0.5, 0},
    {0, 0, 0, 0, 0, 0, 1},
    {0, 0, 0, 0, 0, 0, 1},
    {0, 0, 0, 0, 0, 0, 0}
   } );
fun[t_] := ( {
    {0, f1[1/6, t], 0, 0, 0, 0, 0},
    {0, 0, f2[1/2, t], 0, 0, 0, 0},
    {0, f6[1/1, t], 0, f4[1/1, t], f3[1/1, t], 0, 0},
    {0, f7[1/5, t], 0, 0, f5[1/2, t], f8[1/8, t], 0},
    {0, 0, 0, 0, 0, 0, f9[1, t]},
    {0, 0, 0, 0, 0, 0, f10[1, t]},
    {0, 0, 0, 0, 0, 0, 0}
   } );
q[t_] := p*fun[t];
A = Table[
   Table[q[\[Nu]*h][[k, j]] + 
     Sum[(q[(\[Nu] - 1)*h][[i, j]] - 
          q[((\[Nu] - 1)*h + \[Nu]*h)/2][[i, j]])*
        Subscript[\[Phi], i][(\[Nu] - \[Nu])*h] + 
       Sum[(q[((\[Nu] - \[Mu])*h + (\[Nu] - \[Mu] + 1)*h)/2][[i, j]] -
            q[((\[Nu] - \[Mu] - 1)*h + (\[Nu] - \[Mu])*h)/2][[i, j]])*
         Subscript[\[Phi], i][h*(\[Mu] - 1)], {\[Mu], 1, \[Nu] - 1}] +
        q[h/2][[i, j]]*Subscript[\[Phi], i][h*(\[Nu] - 1)], {i, 2, 
       7}], {j, 1, 7}, {k, 1}], {\[Nu], 1, 100}];

A // MatrixForm
B = Table[A[[\[Nu]]], {\[Nu], 1, 2}];
B // MatrixForm
aa[t_] := 
  Table[Subscript[\[Phi], j][h*\[Nu]] == A[[\[Nu], j]], {\[Nu], 1, 
    3}, {j, 1, 7}];
bb[t_] := 
  Table[Subscript[\[Phi], j][h*\[Nu]], {\[Nu], 1, 3}, {j, 1, 7}];
aa[t] // MatrixForm

bb[t]
sysLT[t_] := 
  Table[Solve[aa[t][[\[Nu]]], bb[t][[\[Nu]]]], {\[Nu], 1, 3}];

sysLT[t][[2]] // MatrixForm

aa[t_] := Table[bb[t][[\[Nu]]] /. sysLT[[\[Nu]]], {\[Nu], 1, 3}];

aa[t] // MatrixForm