error: index exceeds matrix dimensions

I have the following code that I am trying to run but keep getting the "index exceeds matrix dimensions" error for the rho_a line. Does anyone know how to fix this error?

h1=0;r=0;go=9.81;Mi_1=298000;

for t=0:100; 

      Pc2_RD_1(t+1) = Pc_RD*(1+(.015*dt_1)^2);
      Pe_RD_1(t+1) = Pc2_RD_1(t+1)*(1+((gamma_RD-1)/2)*Me_RD^2)^(gamma_RD/(1-gamma_RD));
      Isp_RD_1(t+1) = (e_RD/Pc2_RD_1(t+1))*(Pe_RD_1(t+1)-Pa));
      mprop_RD1(t+1) = At_RD*Pc2_RD_1(t+1)/c_RD;
      thrust_RD1(t+1) = Isp_RD_1(t+1)*g0*mprop_RD1(t+1)/1000;
      dM_RD1(t+1) = mprop_RD1(t+1)*dt_1
      dU1(t+1) = Isp_RD_1(t+1)*g0*dM_RD1(t+1)/Mi_1
      rho_a(t+1) = 1.2*exp(((-2.9e-5)*h1(t+1)^1.15))
      D(t+1) = .5*C_D*Af*rho_a(t+1)*dU1(t+1)^2
      g(t+1) = g0*(r_E/(r_E+h1))^2
      Uy_dt(t+1) = dU1(t+1)*cos(theta)-dt_1(t+1)*(D(t+1)/(Mi_1-dM_RD1(t+1))+g(t+1))*cos(theta)
      Ux_dt(t+1) = dU1(t+1)*sin(theta(t+1))-dt_1(t+1)*(D(t+1)/(Mi_1-dM_RD1(t+1)))*sin(theta(t+1))
      U(t+1) = (Ux_dt(t+1))^2+(Uy_dt(t+1))^2
      h1(t+1) = h1(t+1)+Uy_dt(t+1)*dt_1(t+1)
      r(t+1) = r(t+1)+Ux_dt(t+1)*dt_1(t+1)
      if Ux_dt(t+1)>0
          theta(t+1) = atan(Ux_dt(t+1)/Uy_dt(t+1))
      else
          theta(t+1) = atan(Ux_dt(t+1)/Uy_dt(t+1))+pi()
      end
      dt_1(t+1)=dt_1+.2;
      time(t+1)=t+.2;
  end

[Image Analyst]

That's not the error we get. We get

Error: File: test3.m Line: 16 Column: 63
Unbalanced or unexpected parenthesis or bracket.

Please post your actual code so that we can run it. Also, read this: http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

[the cyclist]

Also, we don't have the definitions of your variables such as

Pc_RD

so, as Image Analyst suggests, seeing your actual code will help us help you.

[Shannon Hemp]

Ru = 8314; % J/(kmol*K)
g0 = 9.81; % m/s^2
C_D = 0.5;
r_E = 6378e3; % m
Pa = 1.01325e5; % Pa
h1 = 0;
U = 0;
r = 0;
theta = 0;
phi = pi()/4;
Mi_1 = 284089; % kg

Df = 4.2; % m
Af = (Df^2)*pi()/4; 

Tc_RD = 3600; % K
M_RD = 23.6; % kg/kmol
gamma_RD = 1.217;
Pc_RD = 25.662e6; % Pa
De_RD = 1.43; % m
eta_RD = 1.056;
lamda_RD = 0.98;
e_RD = 36.4;
At_RD = .08824; % m^2

Me_RD = 4.1;

c_RD = (eta_RD*sqrt(gamma_RD*(Ru/M_RD)*Tc_RD))/(gamma_RD*((2/(gamma_RD+1))^((gamma_RD+1)/(2*(gamma_RD-1)))));

dt_1=0.2;
for t=0:100;
    Pc2_RD_1(t+1) = Pc_RD*(1+(.015*dt_1)^2)
    Pe_RD_1(t+1) = Pc2_RD_1(t+1)*(1+((gamma_RD-1)/2)*Me_RD^2)^(gamma_RD/(1-gamma_RD));
    Isp_RD_1(t+1) = (lamda_RD*c_RD/g0)*(gamma_RD*sqrt((2/(gamma_RD-1))*((2/(gamma_RD+1))^((gamma_RD+1)/(gamma_RD-1)))*(1-((Pe_RD_1(t+1)/Pc2_RD_1(t+1))^((gamma_RD-1)/gamma_RD))))+(e_RD/Pc2_RD_1(t+1))*(Pe_RD_1(t+1)-Pa));
    mprop_RD1(t+1) = At_RD*Pc2_RD_1(t+1)/c_RD;
    thrust_RD1(t+1) = Isp_RD_1(t+1)*g0*mprop_RD1(t+1)/1000;
    dM_RD1(t+1) = mprop_RD1(t+1)*dt_1
    dU1(t+1) = Isp_RD_1(t+1)*g0*dM_RD1(t+1)/Mi_1
    rho_a(t+1) = 1.2*exp(((-2.9e-5)*h1(t+1)^1.15))
    D(t+1) = .5*C_D*Af*rho_a(t+1)*dU1(t+1)^2
    g(t+1) = g0*(r_E/(r_E+h1))^2
    Uy_dt(t+1) = dU1(t+1)*cos(theta)-dt_1(t+1)*(D(t+1)/(Mi_1-dM_RD1(t+1))+g(t+1))*cos(theta)
    Ux_dt(t+1) = dU1(t+1)*sin(theta(t+1))-dt_1(t+1)*(D(t+1)/(Mi_1-dM_RD1(t+1)))*sin(theta(t+1))
    U(t+1) = (Ux_dt(t+1))^2+(Uy_dt(t+1))^2
    h1(t+1) = h1(t+1)+Uy_dt(t+1)*dt_1(t+1)
    r(t+1) = r(t+1)+Ux_dt(t+1)*dt_1(t+1)
    if Ux_dt(t+1)>0
        theta(t+1) = atan(Ux_dt(t+1)/Uy_dt(t+1))
    else
        theta(t+1) = atan(Ux_dt(t+1)/Uy_dt(t+1))+pi()
    end
    dt_1(t+1)=dt_1+.2;
    time(t+1)=t+.2;
end

[Image Analyst]

This uncommented code is impenetrable. It looks like alphabet soup to me. I suggest you step through it a line at a time in the debugger to figure out what it's doing.

Answer 1

[Image Analyst]

You set h1 to be zero, a scalar. So why do you think there should be a (t+1)'th element as if it were an array? It's not an array.

[Shannon Hemp]

Ok I see what you're saying. What my problem is is I initially want to solve the rho at h1=0 and have that go through the iteration. Then at the end I want to find a new h1 and go through the iteration again. For some reason I can't figure out how to do this.

Answer 2

[Guillaume]

What's puzzling is why you can't see the problem.

rho(t+1) depends on h1(t+1), which does not exists yet, hence why you get the error. h1(t+1) is calculated later but since it depends on rho(t+1), you've got a circular reference. Assuming the circular reference is a mistake, only you can tell which of the equations is wrong.

[Shannon Hemp]

Ok I see what you're saying. What my problem is is I initially want to solve the rho at h1=0 and have that go through the iteration. Then at the end I want to find a new h1 and go through the iteration again. For some reason I can't figure out how to do this.

[Guillaume]

The problem is that you're using t both as an index and as your time. Ideally, the two should be decoupled, and you'd have an index going from 1 to numel(t) with t a vector of times (e.g 0:100 for example).

t = 0:100;
h1 = [h0, zeros(size(t))]; %h1(1) is initial condition
PC2_RD = [PC_RD, zeros(size(t))]; %PC_RD is initial condition
%... etc. for all other arrays. They all should have numel(t)+1 elements
for iter = 1:numel(t)
   PC2_RD(iter+1) = PC2_RD(iter) * (1+(.015*dt_1)^2); %is dt_1 supposed to change during the iteration?
   %... 
   rho_a(iter+1) = 1.2*exp(((-2.9e-5)*h1(iter)^1.15))
   %...
end