Why am I getting error message "array indices must be positive integers or logical values"

When I run my code I keep getting the error "array indices must be positive integers or logical values" and I believe the issues stems from my calculation of mu. However, I have looked at the math and mu should be a scalar as it is a function of Y and thetahatdot, Y is a 1x3 and thetahatdot is a 3x1 because it is a function of scalars .* Yt which is a 3x1. My code is posted below. Any help is much appreciated.

low=0;%stated that all given constants m,c,K,g,a,phi were positive, restricted values of xd and x to being positive only as well
high=5;%restricted upper bound values of constants, xd, and x to 5 as it proved to prevent simulation from "blowing up"
Am1=zeros(1,1)+((high-low).*rand(99,1) + low);
Am2=zeros(1);
Am=vertcat(Am2,Am1);
Ac1=zeros(1,1)+((high-low).*rand(99,1) + low);
Ac2=zeros(1);
Ac=vertcat(Am2,Am1);
AK1=zeros(1,1)+((high-low).*rand(99,1) + low);
AK2=zeros(1);
AK=vertcat(AK2,AK1);
Aa=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant a
Aalpha=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant alpha
lowphi=-1*pi;%ensures lower bound of phi will be -pi so that there is only one equilibrium point
highphi=pi;%ensures upper bound of phi will be pi so that there is only one equilibrium point
Aphi=(highphi-lowphi).*rand(100,1) + lowphi;%column vector of 100 uniformly distributed decimals between -pi and pi for constant phi
AB=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant Beta
AJ=(high-low).*rand(100,1) + low;
Agamma=(high-low).*rand(100,1) + low;
Aga=(high-low).*rand(100,1) + low;
Axd=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 5 for variable xd
xa=zeros(1,1)+((high-low).*rand(99,1) + low);
xb=zeros(1);
Ax=vertcat(xb,xa);%column vector of 100 uniformly distributed decimals between 0 and 5 for variable x with the initial condition x(0)=0
ahat1=zeros(1,1)+((high-low).*rand(99,1) + low);
ahat2=zeros(1);
Aahat=vertcat(ahat2,ahat1);
for ii=1:100%signifies 100 iterations for each variable xd and x, and constants c,K,a,alpha,phi, and B
    m=Am(ii);%m takes on a specific value from the 100 randomly generated values for each iteration 
    c=Ac(ii);%c takes on a specific value from the 100 randomly generated values for each iteration
    K=AK(ii);%K takes on a specific value from the 100 randomly generated values for each iteration
    a=Aa(ii);%a takes on a specific value from the 100 randomly generated values for each iteration
    alpha=Aalpha(ii);%alpha takes on a specific value from the 100 randomly generated values for each iteration
    phi=Aphi(ii);%phi takes on a specific value from the 100 randomly generated values for each iteration
    B=AB(ii);%B takes on a specific value from the 100 randomly generated values for each iteration
    J=AJ(ii);
    gamma=Agamma(ii);
    ga=Aga(ii);
    xd=Axd(ii);%xd takes on a specific value from the 100 randomly generated values for each iteration
    x=Ax(ii);%x takes on a specific value from the 100 randomly generated values for each iteration
    thetahat=[m;c;K];
    ahat=Aahat(ii);
    g=9.81;%gravitational constant in m/s^2
    tspan = [0 10];%time span for elapsed simulation
    options = odeset('reltol', 1e-5, 'abstol', 1e-8 );%default options for relative and absolute tolerance.  The tolerances are used to limit the local discretization error during integration
    x0 = [0 0 0 0];%initial conditions for x,xdot,xdotdot and xdotdotdot
    [t, y] = ode15s(@(t,xdot) dynamics(t,xdot,x,xd,m,c,K,a,alpha,phi,B,J,g,gamma,ga,thetahat,ahat), tspan, x0, options );%implementation of ode function in matlab to plot our system dynamics described in the function beneath the end of this for loop
    figure(1)
    plot(t,y(:,1))%plotting of e(t)
    xlabel('Time')
    ylabel('e(t)')
    hold on
    figure(2)
    plot(t,y(:,2))%plotting of r(t)
    xlabel('Time')
    ylabel('r(t)')
    hold on
    figure(3)
    plot(t,y(:,3))%plotting of utilde(t)
    xlabel('Time')
    ylabel('utilde(t)')
    hold on
    figure(4)
    plot(t,y(:,4))%plotting of mu(t)
    xlabel('Time')
    ylabel('mu(t)')
    hold on
    figure(5)
    plot(t,y(:,5))%plotting of thetatilde(t)
    xlabel('Time')
    ylabel('thetatilde(t)')
    hold on
    figure(6)
    plot(t,y(:,6))%plotting of atilde(t)
    xlabel('Time')
    ylabel('atilde(t)')
    hold on
end
hold off
function system = dynamics(t,xdot,x,xd,m,c,K,a,alpha,phi,B,J,g,gamma,ga,thetahat,ahat)
system = zeros(4,1);%initialize system matrix
system(1) = -1*xdot(1);%first row of system matrix defined by equation for edot
system(2) = -1*xdot(2)+alpha*(-1*xdot(1));%second row of system matrix defined by equation for rdot
Y=[g*sin(phi)+alpha*(-1*xdot(1)) xdot(1) x];
Yt=[g*sin(phi)+alpha*(-1*xdot(1));xdot(1);x];
ud=Y*thetahat+(xd-x)+B*(-1*xdot(1)+alpha*(xd-x));
thetahatdot=gamma*(-1*xdot(1)+alpha*(xd-x)).*Yt;
uddot=Y*thetahatdot+(-1*xdot(1))+B*(-1*xdot(2)+alpha*-1*xdot(1));
u=m*xdot(2)+c*xdot(1)+K*x-m*g*sin(phi);
utilde=ud-u;
mu=Y*thetahatdot+(-1*xdot(1))+B(-1*xdot(2)+alpha*(-1*xdot(1)))+ahat*ud-ahat*utilde+J*utilde+(-1*xdot(1)+alpha*(xd-x));
system(3)=uddot+a*u-mu;%utildedot
thetahatdotdot=gamma.*Yt.*(-1*xdot(2)+alpha*-1*xdot(1));
system(4)=Y*thetahatdotdot+(-1*xdot(2))+B*(-1*xdot(3)+alpha*(-1*xdot(2)))+ahat*uddot-ahat*(uddot+a*u-mu)+J*(uddot+a*u-mu)+(-1*xdot(2)+alpha*(-1*xdot(1)));%mudot
system(5)=-1*thetahatdot;%thetatildedot
system(6)=ga*(utilde)^(2)-ga*utilde*ud;%atildedot
end

[DGM]

The error comes from this subexpression:

B(-1*xdot(2)+alpha*(-1*xdot(1)))

Since B is a scalar, I have to ask if this is supposed to be an indexing operation or a multiplication operation?

[Matthew Tortorella]

You are correct that this is the issue. I went through and added a * and it worked. Thanks