ODE Solver Function File Help - Must return a column vector

I keep trying to run an ode23tb but I keep getting a function file error saying my function file must return a column vector

function h=Julia_Prob1a_fun(t,y)
global  T0 cA0 cB0 cC0 Pe ar PeT Le Nu mc mT N M dsi deta A B
k0=700;
Ea=111e3;           %J/mol
vA=1;
vB=0;
vC=-0.36;
vD=-0.09;
nv=(N-2)*(M-2);
cA=y(1:nv);
cB=y(nv+1:2*nv);
cC=y(2*nv+1:3*nv);
T=y(3*nv+1:4*nv);
bA=zeros(nv,nv);
bB=bA;
bC=bA;
bT=bA;

for i=1:nv
    lnK=5693.5/T(i)+1.077*log(T(i))+5.44e-4*T(i)-1.125e-7*T(i)^2-49170/T(i)^2-13.148;
    K=exp(lnK);
    RT=0.082*T(i)*1e-3;
    r=k0*exp(-Ea/(8.314*T(i)))*(cA(i)*RT)^vA*(cB(i)*RT)^vB*(cC(i)*RT)^vC*(cC(i)*RT)^vD*(1-cC(i)*cC(i)/(cA(i)*cB(i)*K));

    bA(i)=-mc*r;
    bB(i)=-mc*r;
    bC(i)=mc*r;
    bT(i)=mT*r;
    
alfa=Pe/(2*dsi)+1/dsi^2;
gama=-Pe/(2*dsi)+1/dsi^2;
alfaT=Le*(PeT/(2*dsi)+1/dsi^2);
gamaT=Le*(-PeT/(2*dsi)+1/dsi^2);
end

for j=1:(M-2)
    eps=ar^2*(1/(2*j*deta^2)+1/deta^2);
    thta=eps*-1;
    epsT=Le*eps;
    thtaT=Le*thta;
    for i=1:(N-2)
        I=(N-2)*(j-1)+i;
        if i==1
            bA(I)=bA(I)+alfa(j)*cA0;
            bB(I)=bB(I)+alfa(j)*cB0;
            bC(I)=bC(I)+alfa(j)*cC0;
            bT(I)=bT(I)+alfaT(j)*T0;
        end
        if i==(N-2)
            bA(I)=bA(I)+gama(j)*cA(I);
            bB(I)=bB(I)+gama(j)*cB(I);
            bC(I)=bC(I)+alfa(j)*cC0;
            bT(I)=bT(I)+alfaT(j)*T0;
        end
        if i==(N-2)
            bA(I)=bA(I)+gama(j)*cA(I);
            bB(I)=bB(I)+gama(j)*cB(I);
            bC(I)=bC(I)+gama(j)*cC(I);
            bT(I)=bT(I)+gamaT(j)*T(I);
        end
        if j==1
            bA(I)=bA(I)+thta*cA(I);
            bB(I)=bB(I)+thta*cB(I);
            bC(I)=bC(I)+thta*cC(I);
            bT(I)=bT(I)+thtaT*T(I);
        end
        if j==(M-2)
            bA(I)=bA(I)+eps*cA(I);
            bB(I)=bB(I)+eps*cB(I);
            bC(I)=bC(I)+eps*cC(I);
            bT(I)=bT(I)+epsT*(T(I)+Nu*deta*T0)/(1+Nu*deta);
        end       
    end
end
hA=A*cA+bA;
hB=A*cB+bB;
hC=A*cC+bC;
hT=B*T+bT;
h=[hA;hB;hC;hT];
end 

Main File

clc;
clear all;
global T0 cA0 cB0 cC0 tf Pe ar PeT Le Nu mc mT N M dsi deta A B
vm=0.066;       %m/s
D=3.19e-7;      %m2/s
P=1.013;        % bar
rhob=21.6e3;    %gm/m3
rhog=0.37e3;    %gm/m3
kT=38.2/3600;   %J/(m.s.K)
Cp=1;           %J/(gm.K)
delH=-41.1e3;   %J/mol;
UT=180/3600;    %J/(m2.s.K)
L=70e-3;        %m
R=6.35e-3;      %m

ra=linspace(0,R);  %various distances from the center of the pipe to R
vz=vm.*(1.-ra.^2./R.^2);

T0=773;                 %K
cT=P/(0.082*T0)*1e3 ;   %mol/m3
cA0=0.27*cT;            %CO;
cB0=0.72*cT;            %H2O;
cC0=0;                  %CO2;
ar=L/R;
Pe=vz.*L/D;
PeT=rhog*Cp*vz.*L/kT;
Le=kT/(rhog*Cp*D);
Nu=UT*R/kT;
mc=L^2*rhob/D;
mT=Le*(-delH)*rhob*L^2/kT;
N=21;                   %z direction
M=21;                   %r direction
dsi=1/(N-1);            %step in z direction
deta=1/(M-1);           %step size in r direction
si=meshgrid(0:dsi:1);   %dimensionless z
eta=meshgrid(0:deta:1)';%dimensionless eta
nv=(N-2)*(M-2);

alfa=Pe/(2*dsi)+1/dsi^2;
beta=-2/dsi^2-2*ar^2/deta^2;
gama=-Pe/(2*dsi)+1/dsi^2;
alfaT=Le*(PeT/(2*dsi)+1/dsi^2);
betaT=Le*(-2/dsi^2-2*ar^2/deta^2);
gamaT=Le*(-PeT/(2*dsi)+1/dsi^2);
A=zeros(nv,nv);
B=A;

for j=1:(M-2)
    eps=ar^2*(1/(2*j*deta^2)+1/deta^2);
    thta=ar^2*(-1/(2*j*deta^2)+1/deta^2);
    epsT=Le*ar^2*(1/(2*j*deta^2)+1/deta^2);
    thtaT=Le*ar^2*(-1/(2*j*deta^2)+1/deta^2);
    for i=1:(N-2)
        I=(N-2)*(j-1)+i;
        A(I,I)=beta;
        B(I,I)=betaT;
        if i>1
            A(I,I-1)=alfa(i);
            B(I,I-1)=alfaT(i);
        end
        if i<(N-2)
            A(I,I+1)=gama(i);
            B(I,I+1)=gamaT(i);
        end
        if j>1
            A(I,I-(N-2))=thta;
            B(I,I-(N-2))=thtaT;
        end
        if j<(M-2)
            A(I,I+(N-2))=eps;
            B(I,I+(N-2))=epsT;
        end       
    end
end

nv=(N-2)*(M-2);
cAi=cA0*ones(nv,1);
cBi=cB0*ones(nv,1); 
cCi=1e-7*ones(nv,1);
Ti=T0*ones(nv,1);
yin=[cAi;cBi;cCi;Ti];

options = odeset('AbsTol',1e-8,'RelTol',1e-5,'MaxStep',0.1);
tf=1e-4;
tau=linspace(0,tf,100);
[t,y]=ode23tb(@Julia_Prob1a_fun,tau,yin',options);

cA=y(:,1:nv);
cB=y(:,nv+1:2*nv);
cC=y(:,2*nv+1:3*nv);
T=y(:,3*nv+1:4*nv);

cAf=zeros(M,N,length(t));
cBf=cAf;
cCf=cAf;
Tf=cAf;

for j=2:M-1
    for i=2:N-1
        I=(N-2)*(j-2)+(i-1);
        for k=1:length(t)
            cAf(j,i,k)=cA(k,I);
            cBf(j,i,k)=cB(k,I);
            cCf(j,i,k)=cC(k,I);
            Tf(j,i,k)=T(k,I);
        end
    end
end

cAf(1,:,:)=cAf(2,:,:);
cAf(M,:,:)=cAf(M-1,:,:);
cAf(:,N,:)=cAf(:,N-1,:);
cAf(:,1,:)=cA0;
cAL=mean(squeeze(cAf(:,N,:)),1);

Xf=100.*(1-cAf./cA0);
XL=mean(squeeze(Xf(:,N,:)),1);

cBf(1,:,:)=cBf(2,:,:);
cBf(M,:,:)=cBf(M-1,:,:);
cBf(:,N,:)=cBf(:,N-1,:);
cBf(:,1,:)=cB0;
cBL=mean(squeeze(cBf(:,N,:)),1);

cCf(1,:,:)=cCf(2,:,:);
cCf(M,:,:)=cCf(M-1,:,:);
cCf(:,N,:)=cCf(:,N-1,:);
cCf(:,1,:)=cC0;
cCL=mean(squeeze(cCf(:,N,:)),1);

Tf(1,:,:)=Tf(2,:,:);
Tf(M,:,:)=(Tf(M-1,:,:)+Nu*deta*T0)./(1+Nu*deta);
Tf(:,N,:)=Tf(:,N-1,:);
Tf(:,1,:)=T0;
TL=mean(squeeze(Tf(:,N,:)),1);

Answer 1

[Cris LaPierre]

Check the dimensions of h, your output from Julia_Prob1a_fun. The ode solver expects a column vector, which has dimensions mx1. I suspect h is a matrix.