How can I solver this error 'Matrix dimensions must agree.' in ode15s with mass matrix?

This is my code

function MMasa_Columna_ELL_Definitivo clc; close all; load('variables') dp1 = polyder(p1); dp3 = polyder(p3); dp2 = polyder(p2); dp4 = polyder(p4); F = 10; x3F = 0; x1F = 0.45; s0 = 3.8; x3s0 = 1; x1s0 = 0; T = 1; %Tiempo de residencia en cada etapa teórica = [min] O = T/60; %Tiempo de residencia en cada etapa teórica = [horas]
x0 = zeros(12,1); %Condiciones iniciales = moles de i/moles totales tspan = [0,1]; %Tiempo de integración = [horas]
Opciones = odeset('MassSingular',@(t,x)Mmasa(t,x,p1,p2,p3,dp1,dp2,dp3),'MStateDependence','strong'); [t,x] = ode23t(@(t,x)modelo(t,x,s0,x3s0,x1s0,F,x3F,x1F,p1,p2,p3,O),tspan,x0,Opciones); plot(t,x(:,:) end function Mdxdt = modelo(~,x,s0,x3s0,x1s0,F,x3F,x1F,p1,p2,p3,O) %Variables de estado x3s1 = x(1); x3s2 = x(2); x3s3 = x(3); x3s4 = x(4); s1 = x(5); s2 = x(6); s3 = x(7); s4 = x(8); r1 = x(9); r2 = x(10); r3 = x(11); r4 = x(12);
%Ecuaciones auxiliares a=zeros(1,length(p1)); a(1,end)=polyval(p3,x3s1); ffx3s1=roots(p1-a); ffx3s1 = ffx3s1(imag(ffx3s1)==0); ffx3s1=ffx3s1(ffx3s1>0 & ffx3s1<0.25);
b=zeros(1,length(p1)); b(1,end)=polyval(p3,x3s2); ffx3s2=roots(p1-b); ffx3s2 = ffx3s2(imag(ffx3s2)==0); ffx3s2=ffx3s2(ffx3s2>0 & ffx3s2<0.25);
c=zeros(1,length(p1)); c(1,end)=polyval(p3,x3s3); ffx3s3=roots(p1-c); ffx3s3 = ffx3s3(imag(ffx3s3)==0); ffx3s3=ffx3s3(ffx3s3>0 & ffx3s3<0.25);
d=zeros(1,length(p1)); d(1,end)=polyval(p3,x3s4); ffx3s4=roots(p1-d); ffx3s4 = ffx3s4(imag(ffx3s4)==0); ffx3s4=ffx3s4(ffx3s4>0 & ffx3s4<0.25);
fsx3s1 = polyval(p2,x3s1); fsx3s2 = polyval(p2,x3s2); fsx3s3 = polyval(p2,x3s3); fsx3s4 = polyval(p2,x3s4);
frx3s1 = polyval(p3,x3s1); frx3s2 = polyval(p3,x3s2); frx3s3 = polyval(p3,x3s3); frx3s4 = polyval(p3,x3s4);
%Ecuaciones diferenciales dn13dt = ((F*x3F + s2*x3s2 - s1*x3s1 - r1*ffx3s1)/O); dn23dt = ((r1*ffx3s1 + s3*x3s3 - s2*x3s2 - r2*ffx3s2)/O); dn33dt = ((r2*ffx3s2 + s4*x3s4 - s3*x3s3 - r3*ffx3s3)/O); dn43dt = ((r3*ffx3s3 + s0*x3s0 - s4*x3s4 - r4*ffx3s4)/O); dn11dt = ((F*x1F + s2*fsx3s2 - s1*fsx3s1 - r1*frx3s1)/O); dn21dt = ((r1*frx3s1 + s3*fsx3s3 - s2*fsx3s2 - r2*frx3s2)/O); dn31dt = ((r2*frx3s2 + s4*fsx3s4 - s3*fsx3s3 - r3*frx3s3)/O); dn41dt = ((r3*frx3s3 + s0*x1s0 - s4*fsx3s4 - r4*frx3s4)/O); dn1dt = ((F + s2 - s1 - r1)/O); dn2dt = ((r1 + s3 - s2 - r2)/O); dn3dt = ((r2 + s4 - s3 - r3)/O); dn4dt = ((r3 + s0 - s4 - r4)/O); Mdxdt = [dn13dt; dn23dt; dn33dt; dn43dt; dn11dt; dn21dt; dn31dt; dn41dt; dn1dt; dn2dt; dn3dt; dn4dt];
end 
function M=Mmasa(~,x,p1,p2,p3,dp1,dp2,dp3) %Variables de estado x3s1 = x(1); x3s2 = x(2); x3s3 = x(3); x3s4 = x(4); s1 = x(5); s2 = x(6); s3 = x(7); s4 = x(8); r1 = x(9); r2 = x(10); r3 = x(11); r4 = x(12); %Ecuaciones auxiliares c1=s1+r1*(polyval(dp3,x3s1)/polyval(dp1,x3s1)); c2=s2+r2*(polyval(dp3,x3s2)/polyval(dp1,x3s2)); c3=s3+r3*(polyval(dp3,x3s3)/polyval(dp1,x3s3)); c4=s4+r4*(polyval(dp3,x3s4)/polyval(dp1,x3s4));
 c5=s1*polyval(dp2,x3s1)+r1*polyval(dp3,x3s1);
 c6=s2*polyval(dp2,x3s2)+r2*polyval(dp3,x3s2);
 c7=s3*polyval(dp2,x3s3)+r3*polyval(dp3,x3s3);
 c8=s4*polyval(dp2,x3s4)+r4*polyval(dp3,x3s4);

 a=zeros(1,length(p1)); a(1,end)=polyval(p3,x3s1);
 ffx3s1=roots(p1-a); ffx3s1 = ffx3s1(imag(ffx3s1)==0);
 ffx3s1=ffx3s1(ffx3s1>0 & ffx3s1<0.25);

 b=zeros(1,length(p1)); b(1,end)=polyval(p3,x3s1);
 ffx3s2=roots(p1-b); ffx3s2 = ffx3s2(imag(ffx3s2)==0);
 ffx3s2=ffx3s2(ffx3s2>0 & ffx3s2<0.25);

 c=zeros(1,length(p1)); c(1,end)=polyval(p3,x3s1);
 ffx3s3=roots(p1-c); ffx3s3 = ffx3s3(imag(ffx3s3)==0);
 ffx3s3=ffx3s3(ffx3s3>0 & ffx3s3<0.25);

 d=zeros(1,length(p1)); d(1,end)=polyval(p3,x3s1);
 ffx3s4=roots(p1-d); ffx3s4 = ffx3s4(imag(ffx3s4)==0);
 ffx3s4=ffx3s4(ffx3s4>0 & ffx3s4<0.25);

 fsx3s1 = polyval(p2,x3s1);
 fsx3s2 = polyval(p2,x3s2);
 fsx3s3 = polyval(p2,x3s3);
 fsx3s4 = polyval(p2,x3s4);

 frx3s1 = polyval(p3,x3s1);
 frx3s2 = polyval(p3,x3s2);
 frx3s3 = polyval(p3,x3s3);
 frx3s4 = polyval(p3,x3s4);
%Definición de la matriz de masa
  M = [c1 0  0  0  x3s1   0      0      0      ffx3s1 0      0      0
       0  c2 0  0  0      x3s2   0      0      0      ffx3s2 0      0
       0  0  c3 0  0      0      x3s3   0      0      0      ffx3s3 0
       0  0  0  c4 0      0      0      x3s4   0      0      0      ffx3s4
       c5 0  0  0  fsx3s1 0      0      0      frx3s1 0      0      0
       0  c6 0  0  0      fsx3s2 0      0      0      frx3s2 0      0
       0  0  c7 0  0      0      fsx3s3 0      0      0      frx3s3 0
       0  0  0  c8 0      0      0      fsx3s4 0      0      0      frx3s4
       0  0  0  0  1      0      0      0      1      0      0      0
       0  0  0  0  0      1      0      0      0      1      0      0
       0  0  0  0  0      0      1      0      0      0      1      0
       0  0  0  0  0      0      0      1      0      0      0      1];
end

This is the error that I get

Matrix dimensions must agree.
Error in ode23t>itsolve (line 905) rhs = h * feval(odeFcn,t,y,odeArgs{:}) - M * z;
Error in ode23t (line 537) itsolve(Mtnew,tnew,ynew,znew,gamma,h,L,U,P,Q,R,odeFcn,odeArgs,rtol,...
Error in MMasa_Columna_ELL_Def (line 15) [t,x] = ode23t(@(t,x)modelo(t,x,s0,x3s0,x1s0,F,x3F,x1F,p1,p2,p3,O,1),tspan,x0,Opciones);

[darova]

I can't read this code. Can you please format it?

a = 2;
b = 3;
c = a + b;

One operation - one line