nSpecies = 5;
n_ax = 25;
n_rad = 10;
C0 = [0;0;0;0;In];
Y0 = repmat(C0,n_ax,n_rad); % size of Y0 is [125 x 10] (2D matrix)
[tSol,YSol] = ode15s(@packBedFun, [0 50], Y0);
%%
function dY = packBedFun(~,Y)
%% ..
Y = reshape(Y,nSpecies + 1,n_ax,n_rad); % you want to make [6 x 25 x 10] (3D matrix)