Algorith doesnt work in program code for Thomas algorithm

%Perkelties metodas
clc
clear
close all

 A=[3 1 0 0;
     -1 4 3 0;
     0 2 4 -1;
     0 0 2 -3];
 b=[2;-2;1;-1];

 xk=transpose([0 0 0 0]);
 k=0;
 zk=A*xk-b;
 pk=A*xk-b;
 skal=dot(zk,zk);
 nmax=1000;
 err=Inf;
 epsilon=0.01;
 fprintf('k\t\t\t\tx\t\t\t\ty\t\t\t\tz\t\t\t\tu\t\t(zk,zk) \n');
 fprintf('%d\t\t%f\t\t%f\t\t%f\t\t%f\t\t%f \n',k,xk,skal);
 while err> epsilon^2 && k<nmax
    rk=A*pk;
    tauk=(dot(zk,pk))/(dot(rk,pk));
    xk1=xk-tauk*zk;
    zk1=zk-tauk*rk; 
    err=dot(zk1,zk1);
    bk=(dot(zk1,zk1))/(dot(zk,zk));
    pk1=zk1+bk*pk;
    k=k+1;
    xk=xk1;
    zk=zk1;
    pk=pk1;
     
    fprintf('%d\t\t%f\t\t%f\t\t%f\t\t%f\t\t%2.12f\n',k,xk,err);
 end
 fprintf('Tiesinių lygčių sistemų sprendinys : \nx=%f\ny=%f\nz=%f\nu=%f \n',xk1);
fprintf('Prireike %1.0f iteraciju\n',k);

I make code right but I get answer (not convergate):

k				x				y				z				u		(zk,zk) 
0		0.000000		0.000000		0.000000		0.000000		10.000000 
1		1.111111		-1.111111		0.555556		-0.555556		18.086419753086
2		1.179409		-1.691647		0.418959		0.605517		28.861061210084
3		1.544930		-1.568929		1.013909		1.029641		7.691686039965
4		1.698675		-1.528515		1.192588		1.428060		0.450952673023
5		1.534880		-1.699673		1.181806		1.410838		0.209363369705
6		1.493959		-1.728117		1.122120		1.515529		0.327367522911
7		1.515388		-1.686118		1.285408		1.573677		1.648896218295
8		1.508312		-1.642814		1.338964		1.644566		0.645185967859
9		1.511725		-1.646326		1.387652		1.691097		0.168569322009
10		1.540909		-1.851758		1.332724		1.674427		1.265727177535
11		4.338046		-1.095935		-7.216512		-11.399516		62049.823583898433
12		4.323182		-1.182441		-7.154933		-11.246880		156873.820168189093
13		4.315913		-1.225556		-7.124544		-11.171170		12736.281492159262
14		4.310903		-1.254144		-7.104096		-11.120710		260.058031879379
15		4.306398		-1.275432		-7.087987		-11.082771		31.832520633989
16		4.301203		-1.292570		-7.073717		-11.052519		6.968346706510
17		4.293950		-1.307501		-7.059607		-11.027835		2.089329789837
18		4.282115		-1.321920		-7.043601		-11.007902		0.852241344763
19		4.257020		-1.338372		-7.019973		-10.993781		0.683739531413

Which are allways same and I think there is a mistake in my code but I cant find it.

Answer 1

[Steven Lord]

If this is a homework assignment, I'm guessing there may be at least one or two examples in your textbook where they walk through the steps of the algorithm, showing the result of each step. If that's the case, try running through that same example in your code using the debugging tools included in MATLAB. Check the result of each step of your code against the results your textbook says you should receive at that step. Once you've found a discrepancy, determine why the two sets of results differ.

If there is no such example in your textbook, either search online for "thomas algorithm worked example" or work through the steps yourself using pencil and paper (or marker and whiteboard, or other writing implement and writing surface.)

If this is not a homework assignment, consider just using the backslash operator \ instead.