Block Tridiagonal Matlab code
 

Hello

I am trying to write a simple code to solve a block tridiagonal matrix with Tomas method
I developed the code same as a simple tridiagonal matrix but I don't know why my code is going wrong
the Matlab code is attached here
%-------------------------------------------------------------
clear all
clc;
y_0 = 0;
h = 40e-3;
y_h = y_0 + h;
dy = 1e-2;
U0 = 40;
nu = 217e-6;
y = y_0:dy:y_h;
Ny = length(y);
dt_B = 0.01;
t_end = 1.08;
t_B = 0:dt_B:t_end;
Nt_B = length(t_B);
d_B = (nu*dt_B/dy);
v_B1 = ones(1,2*Ny);
for j=1:Ny-1
v_B1(2*j)= -dy/2;
end
v_B1(2*Ny)=0;

v_B2 = ones(1,2*Ny-1);v_B2(1) = 0;
for j=2:Ny
v_B2(2*j-1) = -d_B;
end
v_B3 = -1*ones(1,2*Ny-1);v_B3(2*Ny-1)= 1;
for j=2:Ny
v_B3(2*(j-1)) = d_B;
end
v_B4 = zeros(1,2*Ny-2);
for j=2:Ny
v_B4(2*(j-1)) = -dy/2;
end
v_B5 = ones(1,2*Ny-2);
for j=2:Ny
v_B5(2*(j-1)) = 0;
end

v_B = diag(v_B1) + diag(v_B2,1) + diag(v_B3,-1)+diag(v_B4,2)+diag(v_B5,-2);
v_B_new = diag(ones(2*Ny,1));

uv_B = zeros(2*Ny,Nt_B);
%B.C
uv_B(2*Ny,:)= 0;
uv_B(1,:)= U0;
R_B = zeros(2*Ny,1);
R_B(1) = U0;
R_B(2*Ny)=0;
R_B_new = zeros(2*Ny,1);


k=1;
i=1;
A_old(1,1) = v_B(2*i-1,2*i-1);
A_old(1,2) = v_B(2*i-1,2*i);
A_old(2,1) = v_B(2*i,2*i-1);
A_old(2,2) = v_B(2*i,2*i);
%--------------------------------
R_old(1,1) = U0;
R_old(2,1) = 0;
%--------------------------------
R_new = A_old\R_old;
%--------------------------------
R_B_new(2*i-1) = R_new(1,1);
R_B_new(2*i,1) = R_new(2,1);
%--------------------------------
for k = 1:Nt_B-1
for j=2:Ny-1
R_B(2*j-1)= uv_B(2*j-1,k) + uv_B(2*(j-1)-1,k) + d_B*(uv_B(2*j,k) - uv_B(2*(j-1),k));
end
for i=2:Ny
%i = 2;
%get B
B_old(1,1) = v_B(2*(i-1)-1,2*i-1);
B_old(1,2) = v_B(2*(i-1)-1,2*i);
B_old(2,1) = v_B(2*(i-1),2*i-1);
B_old(2,2) = v_B(2*(i-1),2*i);
%get C
C_old(1,1) = v_B(2*i-1,2*(i-1)-1);
C_old(1,2) = v_B(2*i-1,2*(i-1));
C_old(2,1) = v_B(2*i,2*(i-1)-1);
C_old(2,2) = v_B(2*i,2*(i-1));
%get R
R_old(1,1) = uv_B(2*i-1,k) + uv_B(2*(i-1)-1,k) + d_B*(uv_B(2*i,k) - uv_B(2*(i-1),k));
R_old(2,1) = 0;
A_new(1,1) = v_B(2*i-1,2*i-1);
A_new(1,2) = v_B(2*i-1,2*i);
A_new(2,1) = v_B(2*i,2*i-1);
A_new(2,2) = v_B(2*i,2*i);
%----------------------------------
B_new = A_old\B_old;
%----------------------------------
v_B_new(2*(i-1)-1,2*i-1) = B_new(1,1);
v_B_new(2*(i-1)-1,2*i) = B_new(1,2);
v_B_new(2*(i-1),2*i-1) = B_new(2,1);
v_B_new(2*(i-1),2*i) = B_new(2,2);
%----------------------------------
A_old = A_new - C_old*B_new;
R_new = A_old\(R_old - C_old*R_new);
%--------------------------------
R_B_new(2*i-1) = R_new(1,1);
R_B_new(2*i,1) = R_new(2,1);
end
i = Ny;
uv_B(2*i-1,k+1)= R_B_new(2*i-1);
uv_B(2*i,k+1)= R_B_new(2*i);
for i = Ny-1:-1:1
uv_B(2*i-1,k+1)= R_B_new(2*i-1) - v_B_new(2*i-1,2*(i+1)-1)*R_B_new(2*(i+1)-1)-...
v_B_new(2*i-1,2*(i+1))*R_B_new(2*(i+1));
uv_B(2*i,k+1)= R_B_new(2*i) - v_B_new(2*i,2*(i+1)-1)*R_B_new(2*(i+1)-1)-...
v_B_new(2*i,2*(i+1))*R_B_new(2*(i+1));
end
end

I haven't tried it, but you may want to give a look here:

https://it.mathworks.com/matlabcentr...iagonal-solver