How do I solve a PDE which included composition of functions?

I have a PDE which has to be solved iteratively. The equation is,

with initial condition

Here y represents the iteration number and not a power of second derivative. Each y will have its own set of parameters. That is the vectors, matrices and probability generating functions F(s,t)

So I start with y=1 and solve

Next I want to extract the value for 𝐹1(𝑠,𝑡) 𝑠=[0,...,0].

Then I got to iteration y=2 and solve

In the second iteration,s has been replaced with the probability generating function of the previous iteration, so it takes up the pgf which

started in the first iteration with initial condition that was specified then. So the pgf has evolved. Now I want to find at y=2 for t∈[0,1] the value of

when 𝑠=[0,...,0]. How do I solve this numerically on matlab?

Point to remember is that all methods I have tried including ODE45 solver is taking the solution of the previous iteration which is the value when 𝑠=[0,...,0] and using it in the second iteration. This is not what I need actually. I need to use the previous function in the second iteration and solve it as an independent system of matrix differential equation. Would be great if I could get some help on this?

[Thisaakhya]

I have the code which extracts the matrices for the first two iterations. I need to solve the PDEs after this.

%% Step 1: Load Data
opts = detectImportOptions('rates_data.xlsx');
rates_data = readtable('rates_data.xlsx', opts);

% Filter for Congo
le_data = rates_data(strcmp(rates_data.country, 'Congo'), :);

%% Define constants
n = 101;  % Number of age groups
vector1 = ones(n, 1);
alpha = zeros(1, n);
alpha(1) = 1;  % Initial population

%% Initial condition for first iteration
F0 = zeros(n, 1);
F0(1) = 1;
tspan = [0 1]; 

%% Extract Data for Year 1950 (y=1 iteration)
le_data_1950 = le_data(le_data.year == 1950, :);

% Construct Matrices
D1_1950 = diag(le_data_1950.gamma);  
lambda_1950 = le_data_1950.lambda;
gamma_1950 = le_data_1950.gamma;
mu_1950 = le_data_1950.mu;
d_1950 = mu_1950;

% Initialize D0_1950 matrix
D0_1950 = zeros(n, n);
for j = 1:n
    D0_1950(j, j) = -(lambda_1950(j) + gamma_1950(j) + mu_1950(j));
    if j < n
        D0_1950(j, j+1) = lambda_1950(j);
    end
end

% Compute B_1950
B_1950 = kron(alpha, D1_1950);

%% Extract Data for Year 1951 (y=2 iteration)
le_data_1951 = le_data(le_data.year == 1951, :);

% Construct Matrices for 1951
D1_1951 = diag(le_data_1951.gamma);
lambda_1951 = le_data_1951.lambda;
gamma_1951 = le_data_1951.gamma;
mu_1951 = le_data_1951.mu;
d_1951 = mu_1951;

% Initialize D0_1951
D0_1951 = zeros(n, n);
for j = 1:n
    D0_1951(j, j) = -(lambda_1951(j) + gamma_1951(j) + mu_1951(j));
    if j < n
        D0_1951(j, j+1) = lambda_1951(j);
    end
end

% Compute B2_1951
B_1951 = kron(alpha, D1_1951);