Error in using pdepe function

Hello everyone, I want to solve two partial differential equations (in two different co-ordinate systems) simultaneously. But since it not possible to solve them in a single function, I have written them in two separate functions namely 'solidsolve' and 'fluidsolve' and I am executing these two functions in loop over several timespans.

But, while executing this code, I am getting an error message:

Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 0-by-0.

Error in finalcode>solidic (line 73)

sI(1, 1) = T0;

Error in pdepe (line 229)

temp = feval(ic,xmesh(1),varargin{:});

Error in finalcode>solidsolve (line 42)

sol = pdepe(2,@solidge,@solidic,@solidbc,x,t);

Error in finalcode (line 17)

sols = solidsolve(t);

I already tried the normal way to do it. Since it failed, I tried by first declaring a matrix of zeros in every function and then assigning the values in it. This too failed with the above written error message. What must be done to match the dimensions?

Thanks to the community.

clc 

tstart = 0; 
steps = 10;
tend = 0.01;
interval = (tend - tstart) / steps;
tbegin = tstart;
tstop = tbegin + interval;

Tsi = 300; Tffin = 600;
Csi = 0.001; Cffin = 0;

while tbegin < tend
    
    t = linspace(tbegin, tstop, 3);
         
    sols = solidsolve(t);
    Ts = sols(:,:,1);
    Cs = sols(:,:,2);
    
    Tsval = mean(Ts(3, :));
    Tsfin = Tsval; Tfi = Tsval; Tfb = Tsval;
    Csval = mean(Cs(3, :));
    Csfin = Csval; Cfi = Csval; Cfb = Csval;
    
    solf = fluidsolve(t);
    Tf = solf(:,:,1);
    Cf = solf(:,:,2);
    
    Tfval = mean(Tf(3, :));
    Tsi = Tfval; Tffin = Tfval;
    Cfval = mean(Cf(3, :));
    Csi = Cfval; Cffin = Cfval;
    
    tbegin = tbegin + interval;
    tstop = tstop + interval;
    
end

%solver function for solid
function sol = solidsolve(t)
    x = linspace(0,0.00001,3);
    sol = pdepe(2,@solidge,@solidic,@solidbc,x,t);
end

%governing equation function for solid
function sG = solidge(x,t,u,DuDx)
    sG = zeros(3, 2);
    rho = 3665.26358;
    Cp = 2129.596457;
    Hr = -2341.975;
    k = -391.781;
    %R = k*u(2);
    kf = 0.467442272;
    De = 0.000000002;
    %sG(1) = [rho*Cp/kf; 1/De]; 
    %sG(2) = [1; 1].*DuDx;
    %sG(3) = [-Hr*k*u(2)/kf; k*u(2)/De];
    
    sG(1, 1) = rho*Cp/kf; 
    sG(2, 1) = DuDx(1);
    sG(3, 1) = -Hr*k*u(2)/kf;
    sG(1, 2) = 1/De; 
    sG(2, 2) = DuDx(2);
    sG(3, 2) = k*u(2)/De;
end

%initial condition function for solid
function sI = solidic(x)
    sI = zeros(2, 1);
    global Tsi Csi
    T0 = Tsi;
    C0 = Csi;
    %sI(1) = [T0; C0]; 
    
    sI(1, 1) = T0;
    sI(2, 1) = C0;
end

%boundary condition function for solid
function sB = solidbc(xl,ul,xr,ur,t)
    sB = zeros(4, 2);
    global Tffin Cffin
    h = 1500;
    ks = 0.467442272;
    km = 0.000000002;
    kd = 0.000000002;
    %sB(1) = [0; 0]; 
    %sB(2) = [1; 1]; 
    %sB(3) = [h*(ur(1)-Tffin)/ks; km*(ur(2)-Cffin)/kd]; 
    %sB(4) = [1; 1];
    
    sB(1,1) = 0; 
    sB(2,1) = 1; 
    sB(3,1) = h*(ur(1)-Tffin)/ks; 
    sB(4,1) = 1;
    sB(1,2) = 0; 
    sB(2,2) = 1; 
    sB(3,2) = km*(ur(2)-Cffin)/kd; 
    sB(4,2) = 1;
end

%solver function for fluid
function sol = fluidsolve(t)
    y = linspace(0,2,3);
    sol = pdepe(0,@fluidge,@fluidic,@fluidbc,y,t);
end

%governing equation function for fluid
function fG = fluidge(x,t,u,DuDx)
    fG = zeros(3, 2);
    rho = 1.165;
    Cp = 29.124;
    v = 0.1;
    hs = 50;
    as = 4.18*pow(10, -15);
    hf = 50;
    af = 0.1414;
    km = 0.32;
    Tfurn = 700;
    global Tsfin Csfin
    %fG(1) = [-1/v; -1/v]; 
    %fG(2) = [1; 1].*u;
    %fG(3) = [-1*((hs*as*(u(1)-Tsfin)/(rho*Cp*v)) + (hf*af*(u(1)-Tfurn)/(rho*Cp*v))); -1*km*as*(u(2)-Csfin)];
    
    fG(1, 1) = -1/v; 
    fG(2, 1) = u(1);
    fG(3, 1) = -1*((hs*as*(u(1)-Tsfin)/(rho*Cp*v)) + (hf*af*(u(1)-Tfurn)/(rho*Cp*v)));
    fG(1, 2) = -1/v; 
    fG(2, 2) = u(2);
    fG(3, 2) = -1*km*as*(u(2)-Csfin);
end

%initial condition function for fluid
function fI = fluidic(x);
    fI = zeros(1, 2);
    global Tfi Cfi
    T0 = Tfi;
    C0 = Cfi;
    %fI(1) = [T0; C0];
    
    fI(1, 1) = T0;
    fI(1, 2) = C0;
end

%boundary condition function for fluid
function fB = fluidbc(xl,ul,xr,ur,t)
    fB = zeros(4, 2);
    global Tfb Cfb
    %fB(1) = [Tfb; Cfb]; 
    %fB(2) = [0; 0]; 
    %fB(3) = [0; 0]; 
    %fB(4) = [0; 0];
    
    fB(1, 1) = Tfb; 
    fB(2, 1) = 0; 
    fB(3, 1) = 0; 
    fB(4, 1) = 0;
    fB(1, 2) = Cfb; 
    fB(2, 2) = 0; 
    fB(3, 2) = 0; 
    fB(4, 2) = 0;
end