How can I change this code to form a continuous plotted line?

All I get are invisible plots unless I change then, for example, to circles then I can see them. But this is not very neat nor what I need. What do I need to change?

b_v = 3.7;                     % volts of battery
b_a = 24000e-03;               % 6000 mA x 4 cells
b_r = 3.5e-03;                 % resistance of battery, 1-7 mOhms

t_sun = 1;                     % time in the sun
t_ecl = 1;                     % time in eclipse

t_on  = 800;                   % longest time for battery to be on, s
t_off = 228;                   % shortest time for battery to be on, s

m_2U  = 2.66;                  % mass of a 2U cubesat, kg
shc_b = 879.2;                 % specific heat of battery
shc_g = 130;                   % specific heat capacity of gold

A_s   = 0.1*0.1*2+0.1*0.227*4; % surface area of 2U
p_l   = 1/4;                   % fraction of long side are louvers, 0.1*0.227 sides
A_l   = p_l*0.1*0.227;         % surface of area of louvers
A_sl  = A_s - A_l;             % surface area that are not the louvers

sigma = 5.67e-08;              % Boltzmann constant
e_g   = 0.04;                  % emissivity of gold
e_b   = 1;                     % ideal black body emissivity 

T_cub_n = 333;                 % cubesat starting temperature, kelvin, 50 C
T_cub_l = 333;                 % cubesat starting temperature with louvers, 50 C
%% Variables


for t = 1:20
    
    % without louvers open
    Q_n(t)   = sigma*e_g*A_s*-(T_cub_n(t)^4); % Heat without louvers
    dT_n(t)  = Q_n(t)/(shc_g*m_2U);    % change in temperature per second
    T_cub_n(t+1) = T_cub_n(t) + dT_n(t);
    
    % with louvers open
    Q_l(t)   = sigma*e_b*A_l*-(T_cub_l(t)^4)+sigma*e_g*A_sl*-(T_cub_l(t)^4); %heat with louvers
    dT_l(t)  = Q_l(t)/(shc_g*m_2U);    % change in temperature per second
    T_cub_l(t+1) = [T_cub_l(t) + dT_l(t)];
    
    time = t - 1;
    
    plot(time, T_cub_n(t), 'r')
    xlabel('Time (s)')
    ylabel('Temperature (K)')
    title('Temperature of CubeSat')
    
    hold on 
    plot(time, T_cub_l(t), 'b')

end

Answer 1

[KALYAN ACHARJYA]

Modified Done:

1. Considering time as new variable, t form 1 to 20-1, as t=1 in defined before loop, may be considering initial values.
2. Time is defined just before the loop, as it has no role within the loop, except to get the number of interations.
3. Preallocation of T_cub_n and T_cub_l (Always recomended)
4. No need of array indexing of Q_n, dT_n, Q_l and dT_l, as all those variable just used once to calculate T_cub_n,T_cub_l respectively.

b_v = 3.7;                     % volts of battery
b_a = 24000e-03;               % 6000 mA x 4 cells
b_r = 3.5e-03;                 % resistance of battery, 1-7 mOhms
t_sun = 1;                     % time in the sun
t_ecl = 1;                     % time in eclipse
t_on  = 800;                   % longest time for battery to be on, s
t_off = 228;                   % shortest time for battery to be on, s
m_2U  = 2.66;                  % mass of a 2U cubesat, kg
shc_b = 879.2;                 % specific heat of battery
shc_g = 130;                   % specific heat capacity of gold
A_s   = 0.1*0.1*2+0.1*0.227*4; % surface area of 2U
p_l   = 1/4;                   % fraction of long side are louvers, 0.1*0.227 sides
A_l   = p_l*0.1*0.227;         % surface of area of louvers
A_sl  = A_s - A_l;             % surface area that are not the louvers
sigma = 5.67e-08;              % Boltzmann constant
e_g   = 0.04;                  % emissivity of gold
e_b   = 1;                     % ideal black body emissivity 
time=1:20;
T_cub_n=zeros(1,length(time));
T_cub_l=zeros(1,length(time));
T_cub_n(1)= 333;                 % cubesat starting temperature, kelvin, 50 C
T_cub_l(1)= 333;                 % cubesat starting temperature with louvers, 50 C
%% Variables
for t=1:length(time)-1
    % without louvers open
    Q_n= sigma*e_g*A_s*-(T_cub_n(t)^4); % Heat without louvers
    dT_n= Q_n/(shc_g*m_2U);    % change in temperature per second
    T_cub_n(t+1)=T_cub_n(t) + dT_n;
    
    % with louvers open
    Q_l= sigma*e_b*A_l*-(T_cub_l(t)^4)+sigma*e_g*A_sl*-(T_cub_l(t)^4); %heat with louvers
    dT_l= Q_l/(shc_g*m_2U);    % change in temperature per second
    T_cub_l(t+1)=T_cub_l(t) + dT_l;
end
plot(time, T_cub_n, 'r',time,T_cub_l, 'b')
xlabel('Time (s)')
ylabel('Temperature (K)')
title('Temperature of CubeSat')
grid on;

[Jonathan Charman]

Thank you very much.