HELP! Index in position 3 exceeds array bounds. Index must not exceed 1.

    % Alternating-Direction Implicit Solution Method
    % u[xx]+u[yy]=u[t]
    % Boundary conditions  u(x,0,t) = 0  u(x,1,t) = 1
    %                      u(0,y,t) = 0  u(1,y,t) = 1
    % initial conditions   u(x,y,0) = 0  0<=x<1    0<=y<1
    % i = spatial index in x-direction, from 1 to R
    % j = spatial index in y-direction, from 1 to S

R=21; % last x-point
S=21; % last y-point
N=20; % last time step = N+1

    %% constatnts
dx = 1/(R-1);
dx2 = dx*dx;
dy = 1/(S-1);
dy2 = dy*dy;
dxdy = dx2/dy2;
dydx = dy2/dx2;
dt = dx2; % dt to dx2 for good stability and results
    % independent space variables
x=0:dx:1;
y=0:dy:1;
    %% sizing matrices
u=zeros(R,S,N);
u1=zeros(R,S);
t=zeros(1,N+1);
a=zeros(1,R); 
b=zeros(1,R);
c=zeros(1,R);
d=zeros(1,R);
ba=zeros(1,R);
ga=zeros(1,R);
   %% boundary conditions
for n=1:N
    for i=1:R
        u(i,S,n)=1;
    end
    for j=1:S
        u(R,j,n)=1;
    end
end
    %% intermediate values
for i=1:R
    u1(i,S)=1;
end
for j=1:S
    u1(R,j)=1;
end
    % plot initial conditions
mesh(x,y,u(:,:,1))
xlabel('x-coordinate'); 
ylabel('y-coordinate');
zlabel('u-Temperature');

    %% time step loop
t(1)=0;
for n=1:N
    t(n+1)=t(n)+2*dt;
    % first pass in x-direction
    % first time step - intermediate values at u1(i,j) are calculated
    % constants
dx2dt=dx2/dt;
    % coefficients
for j=2:S-1
    b(2)=-2-dx2dt;
    c(2)=1;
    d(2)=-dxdy*u(2,j-1,n)+(2*dxdy-dx2dt)*u(2,j,n)-dxdy*u(2,j+1,n);
    for i=3:R-2
        a(i)=1;
        b(i)=-2-dx2dt;
        c(i)=1;
        d(i)=-dxdy*u(i,j-1,n)+(2*dxdy-dx2dt)*u(i,j,n)-dxdy*u(i,j+1,n);
    end
    a(R-1)=1;
    b(R-1)=-2-dx2dt;
    d(R-1)=-1-dxdy*u(i,j-1,n)+(2*dxdy-dx2dt)*u(i,j,n)-dxdy*u(i,j+1,n);
        % solution by Thomas Algorithm
    ba(2)=b(2);
    ga(2)=d(2)/b(2);
    for i=3:R-1
        ba(i)=b(i)-a(i)*c(i-1)/ba(i-1);
        ga(i)=(d(i)-a(i)*ga(i-1))/ba(i);
    end
        % back substitution step
    u1(R-1,j)=ga(R-1);
    for i=R-2:-1:2
        u1(i,j)=ga(i)-c(i)*u1(i+1,j,n)/ba(i);
    end
end
    %% second pass in y-direction
    % second time step - final values at u(i,j,n+1) are calculated
    % constants
dy2dt=dy2/dt;
    % coefficients
for i=2:R-1
    b(2)=-2-dy2dt;
    c(2)=1;
    d(2)=-dydx*u1(i-1,2)+(2*dydx-dy2dt)*u1(i,2)-dydx*u1(i+1,2);
    for j=3:S-2
        a(j)=1;
        b(j)=-2-dy2dt;
        c(j)=1;
        d(j)=-dydx*u1(i-1,j)+(2*dydx-dy2dt)*u1(i,j)-dydx*u1(i+1,j);
    end
    a(S-1)=1;
    b(S-1)=-2-dy2dt;
    d(S-1)=-1-dydx*u1(i-1,S-1)+(2*dydx-dy2dt)*u1(i,S-1)-dydx*u1(i+1,S-1);
        % solution by Thomas Algorithm
    ba(2)=b(2);
    ga(2)=d(2)/b(2);
    for j=3:S-1
        ba(j)=b(j)-a(j)*c(j-1)/ba(j-1);
        ga(j)=(d(j)-a(j)*ga(j-1))/ba(j);
    end
        % back substitution step
    u(i,S-1,n+1)=ga(S-1);
    for j=S-2:-1:2
        u(i,j,n+1)=ga(j)-c(j)*u(i,j+1,n+1)/ba(j);
    end
end
end

Index in position 3 exceeds array bounds. Index must not exceed 1.

    %% plot results
mesh(x,y,u(:,:,10))
xlabel('x-coordinate'); 
ylabel('y-coordinate');
zlabel('u-Temperature');

t(10)

mesh(x,y,u(:,:,20))
xlabel('x-coordinate'); 
ylabel('y-coordinate');
zlabel('u-Temperature');
t(20)

Index in position 3 exceeds array bounds. Index must not exceed 1.

Error in deneme1 (line 92)

u1(i,j)=ga(i)-c(i)*u1(i+1,j,n)/ba(i);

what could be the reasaon about this error? I could'n get it what it is. Can you solve and explaain please?

[Jeffrey Clark]

Matrix u1 cannot be referenced as u1(i+1,j,n) which requires it to be of type Multidimensional Arrays - MATLAB & Simulink (mathworks.com)

[Muzaffer Cagislar]

oh okay got it now thanks a lot