Solve System of PDEs with Distributed Initial Conditions

I'm solving a PDE system with initial conditions that are Population Distributed vectored-valued functions for each 6 Populations in the system. I've been using (https://www.mathworks.com/help/matlab/math/partial-differential-equations.html and https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html) as a resource to solve the PDE system with the pdepe() function.

Let's consider the same PDE system from the link provided above, but the Initial condition is adjusted slightly.

Initial Condition:

where and are distrubtion vectored-valued functions that should provide different concentrations of the solutions in different parts of space.

Boundary Condition:

% Discretizing Space and Time:
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1];
t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];

% Solving the PDE System w/ pdepe:
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);

% Extracting Solutions:
u1 = sol(:,:,1);
u2 = sol(:,:,2);

% Plotting the Solution:
figure(1);
surf(x,t,u1)
title('u_1(x,t)')
xlabel('Distance x')
ylabel('Time t')
hold off;

figure(2);
surf(x,t,u2)
title('u_2(x,t)')
xlabel('Distance x')
ylabel('Time t')
hold off;

What adjustments would be made to the Initial Condition function?

Let's establish and to be the following vectored-values distributions between 0 and 1.

delta_1 = rand(13,1);
delta_2 = rand(13,1);

The original code from the link above is

function [c,f,s] = pdefun(x,t,u,dudx) % Equation to solve
c = [1; 1];
f = [0.024; 0.17] .* dudx;
y = u(1) - u(2);
F = exp(5.73*y)-exp(-11.47*y);
s = [-F; F];
end
% ---------------------------------------------
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions
pl = [0; ul(2)];
ql = [1; 0];
pr = [ur(1)-1; 0];
qr = [0; 1];
end
% ---------------------------------------------

But the revisement to the initial condition by setting the initial solutions as distributions is:

function u0 = pdeic(x) % Initial Conditions
n = length(x);
delta_1 = rand(n,1); % u_1(x,0) = delta_1(x)
delta_2 = rand(n,1); % u_2(x,0) = delta_2(x)

u0 = [delta_1'; delta_2'];
end

I got the following results. Are these correct? Is my approach correct?

[Torsten]

pdepe cannot handle unsteady or even stochastic inputs. You will have to use a different solver - whatever your aim might be.

Maybe of interest:

https://en.wikipedia.org/wiki/Stochastic_partial_differential_equation