Error using ==> mrdivide: Matrix dimensions must agree?

[EDIT: 20110609 14:53 CDT - reformat - WDR]

I've been working on an ODE modeling the slipping nature of sand in a cylindrical shell as it rolls down a hill. When I run my code, this error appears:

Error using ==> mrdivide
Matrix dimensions must agree.

Error in ==> sr>f2 at 135
II= (m*a*b*sin(y(1)+alpha)-Iw)/Q+a/J;

Can someone explain to me why this is happening? My initial conditions for this ODE are in a column vector (eg y0 = [PHI;0;0;0];) and I am using ode45 to solve this equation. y0 has been transformed so that it is not confused with a vector (yout = y0.';). Below is the code(note that the sticking system has no errors):

function sr
global T b My AREA a c h STEEL SAND M m alpha Ic Iw g

tstart = 0;
tfinal = 15;

T=0.05; %thickness of wedge
b=0.10; %inner radius
My=2/3*sqrt(T^3*(2*b-T)^3);
AREA= pi*b^2/2-b^2*asin(1-T/b)+(T-b)*sqrt(T*(2*b-T));
a=My/AREA; %centroid calculator (exact)
c=0.11; %outer radius
h=0.55; %height of cylinder

STEEL=7850;
SAND =1602;
M=pi*(c^2-b^2)*h*STEEL;
m=AREA*h*SAND; 
alpha=pi/180*5;
Ic=M/2*(c^2+b^2);
Iw=m*a^2;
g=9.81;
mus=0.5; mud=mus;

PHI=acos((1+M/m)*b/a*sin(alpha));

y0 = [pi/5;0;0;0]; %initial internal angle of wedge, internal angular
%velocity of wedge, position, and velocity of cylinder respectively.
tout = tstart;
yout = y0.'; % use .' at the end so it's not confused with a vector.
teout = [];
yeout = [];
ieout = [];

Q=m*a*sin(y0(1)+alpha)-Ic/b-b*(M+m);

W=Q+m*a*sin(y0(1)+alpha)-Iw/b;
A=1/W*(m*a*g*cos(y0(1)) -(m+M)*b*g*sin(alpha) +m*a*b*y0(2)^2);
PSIDOT=-A/b;

FN=(-A*sin(y0(1)+alpha) - a*PSIDOT + g*cos(y0(1)))/...
    (A*cos(y0(1)+alpha) + a*y0(2)^2 + g*sin(y0(1)));

Z= mus-FN;

for i = 1:30 

options = odeset('Events',@events2,'RelTol',1e-8,'AbsTol',1e-12);
    if tstart == tfinal
        break;
    end;
  [t,y,te,ye,ie] = ode45(@f2,[tstart tfinal],y0,options);

    nt = length(t);
    tout = [tout; t(2:nt)];
    yout = [yout; y(2:nt,:)];

    teout = [teout; te];
    yeout = [yeout; ye];
    ieout = [ieout; ie];

    y0(1) = y(end,1);
    y0(2) = y(end,2);
    y0(3) = y(end,3);
    y0(4) = y(end,4);      
    tstart = t(nt);

options = odeset('Events',@events1,'RelTol',1e-8,'AbsTol',1e-12);
    if tstart == tfinal
        break;
    end;

    [t,y,te,ye,ie] = ode45(@f1,[tstart tfinal],y0,options);

    nt = length(t);
    tout = [tout; t(2:nt)];
    yout = [yout; y(2:nt,:)];

    teout = [teout; te];     
    yeout = [yeout; ye];
    ieout = [ieout; ie];

    y0(1) = y(end,1);
    y0(2) = y(end,2);
    y0(3) = y(end,3);
    y0(4) = y(end,4);
    tstart = t(nt);
  end;

      figure(2)
      subplot(211)
      plot(tout,yout(:,1)) %,teout,yeout(:,1),'ro')
      xlabel('time');
      ylabel('\phi');
      subplot(212)
      plot(tout,yout(:,2)+yout(:,4)/b) %,teout,yeout(:,2),'ro')
      xlabel('time');
      ylabel('d\phi/dt');

      figure(1)
      subplot(211)
      plot(tout,yout(:,3)) %,teout,yeout(:,3),'ro')
      xlabel('time');
      ylabel('x');
      subplot(212)
      plot(tout,yout(:,4)) %,teout,yeout(:,4),'ro')
      xlabel('time');
      ylabel('v');

%for all systems, masses are in kg, lengths are in m, moments of
%intertia are in kg*m^2, and densities are in kg/m^3

function dydt = f1(t,y)
%sticking system
global b a M m alpha Ic Iw g

Q=m*a*sin(y(1)+alpha)-Ic/b-b*(M+m);

W=Q+m*a*sin(y(1)+alpha)-Iw/b;
A=1/W*(m*a*g*cos(y(1)) -(m+M)*b*g*sin(alpha) +m*a*b*y(2)^2*cos(y(1)+alpha));
PSIDOT=-A/b;

dydt = [-y(4)/b;PSIDOT;y(4);A];

function dydt = f2(t,y)
%slipping system
global b a M m alpha Ic Iw g mud

Q=m*a*sin(y(1)+alpha)-Ic/b-b*(M+m);
J=sin(y(1)+alpha)+mud*cos(y(1)+alpha)*sign(y(2)+y(4)/b);

II= (m*a*b*sin(y(1)+alpha)-Iw)/Q+a/J;
PSIDOT=1/II*((g*cos(y(1))-mud*sign(y(2)+y(4)/b)*(g*sin(y(1))+a*y(2)^2))/J-(m*a*g*cos(y(1))-(m+M)*b*g*sin(alpha)+m*a*b*y(2)^2*cos(y(1)+alpha))/Q);

A=1/Q*(m*a*g*cos(y(1))-(m+M)*b*g*sin(alpha)+m*a*b*y(2)^2*cos(y(1)+alpha)+PSIDOT*(m*a*b*sin(y(1)+alpha)-Iw));

dydt = [y(2);PSIDOT;y(4);A];

% --------------------------------------------------------------------------

function [value,isterminal,direction] = events1(t,y)
global b a M m alpha Ic Iw g mus

Q=m*a*sin(y(1)+alpha)-Ic/b-b*(M+m);

W=Q+m*a*sin(y(1)+alpha)-Iw/b;
A=1/W*(m*a*g*cos(y(1)) -(m+M)*b*g*sin(alpha) +m*a*b*y(2)^2);
PSIDOT=-A/b;

FN=(-A*sin(y(1)+alpha) - a*PSIDOT + g*cos(y(1)))/...
    (A*cos(y(1)+alpha) + a*y(2)^2 + g*sin(y(1)));

value = (mus-FN); %finds critical zero
isterminal = 1;   % stops the integration
direction = 0;  

function [value,isterminal,direction] = events2(t,y)
global b
value =y(2)+y(4)/b; 
isterminal = 1;   % stops the integration
direction = 0;

Answer 1

[Walter Roberson]

Change your line

II= (m*a*b*sin(y(1)+alpha)-Iw)/Q+a/J;

to

II= (m.*.a.*.b.*sin(y(1)+alpha)-Iw)./Q + a./J;

I did not trace far enough to see which of your variables is triggering the problem, but if that change gets rid of the division problem then you could put in a breakpoint and display the sizes of the variables.