Index in position 2 is invalid. Array indices must be positive integers or logical values.

%% Material parameters
data_rho = [7840,0.44]; 
data_Cp = [490,0.0733333];   
data_k = [13.1947,0.0126919];     

%% Geometry
L = 0.03;  

%% Boundary conditions
T0 = 1000;      
Tinfinity = 20 + 273;

htc_top = 50;
htc_sides = 100;
htc_btm = 20;
    
Emissivity  = 0.0;   

time = 0;

%% thermocouple positions

tpx = [0.015	0.027	0.027];
tpy = [0.015	0.015	0.027];

%% Simulation control
n = 3; 
CFL = 0.5;

itPlot = 10;

%% Discretize space
dx = L/(n-1);
x = 0:dx:L;
y = L:-dx:0; % T(1,1) is the top left


%% Initialize variables

T = ones(n,n)*T0;
T_new = zeros(n,n);
dtMat = zeros(n,n);
it = 0;

%% Physical constants
stefan = 5.670367e-8; 


%% Record predictions

Temp_thermocouples = interp2(x,y,T,tpx,tpy);

SaveNumber = 1;
TimeVector(SaveNumber) = time;
TempSavedMatrix(SaveNumber,:) = Temp_thermocouples;

%% Model calculation

while max(max(T)) > 300 + 273.15
    %% Thermophysical properties
    rho = data_rho(1) + data_rho(2)*T;
    k0 = data_k(1) + data_k(2)*T;
    Cp = data_Cp(1) + data_Cp(2)*T;
    
    %% Calculate temperature dependent parameters
    
    alpha = k0./rho./Cp;
    Br = stefan*Emissivity*dx./k0;
    Biot_top = htc_top*dx./k0;
    Biot_bottom = htc_btm*dx./k0;
    Biot_sides = htc_sides*dx./k0;
    
    %% time step
    % interior
    for i = 2: n-1
        for j = 2: n-1
            dtMat(i,j) = CFL*dx*dx/4./alpha(i,j);
        end
    end
    
    % top left
    i = 1;
    j = 1;
    dtRadTerm = 0;
    dtTerm = 2+ Biot_top(i,j) + Biot_sides(i,j)+ dtRadTerm;
    dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
    
   
    for j = 2: n-1
         % top edge
        i = 1;
        dtRadTerm = 0;
        dtTerm = 2+ Biot_top(i,j)+ dtRadTerm;
        dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
    end
    
    % top right
    i = 1;
    j = n;
    dtRadTerm = 0;
    dtTerm = 2+ Biot_top(i,j) + Biot_sides(i,j)+ dtRadTerm;
    dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
   
    for i = 2: n-1
       % left edge
        
        j = 1;
        dtRadTerm = 0;
        dtTerm = 2+ Biot_sides(i,j)+ dtRadTerm;
        dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
    
    end
   
    
   
    for i = 2: n-1
         % right edge
        j = n;
        dtRadTerm = 0;
        dtTerm = 2+ Biot_sides(i,j)+ dtRadTerm;
        dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
    
    end
    % btm left
    i = n;
    j = 1;
    dtRadTerm = 0;
    dtTerm = 2+ Biot_bottom(i,j) + Biot_sides(i,j)+ dtRadTerm;
    dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
    
    % btm edge
    for j = 2: n-1
    i = n;
    dtRadTerm = 0;
    dtTerm = 2+ Biot_bottom(i,j)+ dtRadTerm;
    dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
    end
    % btm right
    i = n;
    j = n;
    dtRadTerm = 0;
    dtTerm = 2+ Biot_bottom(i,j) + Biot_sides(i,j)+ dtRadTerm;
    dtMat(i,j) = CFL*dx*dx/2./(alpha(i,j))*dtTerm;
    
    dt = min(dtMat(:));
    
    %% Fourier number
    Fo = alpha*dt/(dx*dx);
    
    %% Calculate the new T vector T^(p+1)
    
    % top left
    
    i = 1;
    j = 1;
    part1 = (1-4*Fo(i,j)*T(i,j));
    part2 = 2*Fo(i,j)*T(i+1,j) + T(i,j+1);
    part3 = 2*Fo(i,j)*(Biot_top(i,j)+Biot_sides(i,j))*(Tinfinity-T(i,j));
    
    T_new(i,j) = part1+part2+part3;
    
   
    for j = 2: n-1
      
         % top surface
        i = 1;
        part1 = (1-4*Fo(i,j)*T(i,j));
        part2a = 2*Fo(i,j)*T(i+1,j);
        part2b = Fo(i,j)*T(i,j+1);
        part2c = Fo(i,j)*T(i,j-1);
        part3 = 2*Fo(i,j)*(Biot_top(i,j))*(Tinfinity-T(i,j));
    
        T_new(i,j) = part1+part2a+part2b+part2c+part3;
    end
    
    % top right
    i = 1;
    j = n;
    part1 = (1-4*Fo(i,j)*T(i,j));
    part2 = 2*Fo(i,j)*T(i+1,j) + T(i,j-1);
    part3 = 2*Fo(i,j)*(Biot_top(i,j)+Biot_sides(i,j))*(Tinfinity-T(i,j));
    
    T_new(i,j) = part1+part2+part3;
    
     for i = 2: n-1
      

         % left surface
        j = 1;
        part1 = (1-4*Fo(i,j)*T(i,j));
        part2a = 2*Fo(i,j)*T(i,j+1);
        part2b = Fo(i,j)*T(i+1,j);
        part2c = Fo(i,j)*T(i-1,j);
        part3 = 2*Fo(i,j)*(Biot_sides(i,j))*(Tinfinity-T(i,j));
    
        T_new(i,j) = part1+part2a+part2b+part2c+part3;
    end
    
     
      
         % right surface
        j = n;
        part1 = (1-4*Fo(i,j)*T(i,j));
        part2a = 2*Fo(i,j)*T(i,j-1);
        part2b = Fo(i,j)*T(i+1,j);
        part2c = Fo(i,j)*T(i-1,j);
        part3 = 2*Fo(i,j)*(Biot_sides(i,j))*(Tinfinity-T(i,j));
    
        T_new(i,j) = part1+part2a+part2b+part2c+part3;
    
    
    % bottom left
    i = n;
    j = 1;
    part1 = (1-4*Fo(i,j)*T(i,j));
    part2 = 2*Fo(i,j)*T(i-1,j) + T(i,j+1);
    part3 = 2*Fo(i,j)*(Biot_bottom(i,j)+Biot_sides(i,j))*(Tinfinity-T(i,j));
    
    T_new(i,j) = part1+part2+part3;
    
    % bottom edge
     %
        i = n;
        part1 = (1-4*Fo(i,j)*T(i,j));
        part2a = 2*Fo(i,j)*T(i-1,j);
        part2b = Fo(i,j)*T(i,j+1);
        part2c = Fo(i,j)*T(i,j-1);
        part3 = 2*Fo(i,j)*(Biot_bottom(i,j))*(Tinfinity-T(i,j));
    
        T_new(i,j) = part1+part2a+part2b+part2c+part3;
   issue is with 
        part2c =
        Fo(i,j)*T(i,j-1);
    

[Ive J]

You are setting j at some point to 1, so j-1 would be zero.

[Jan]

@John McCann: Please post the complete error message and only the relevcant part of the code. Then the readers do not have to guess, where the problem occurs.

Answer 1

[Jan]

The debugger helps you to find the problem. Type in the command window:

dbstop if error

Then run your code again. Matlab will stop when the error occurs and you can examine the values of the locally used variables.