A (and consequently, L, D and U) is a 16x16 matrix. b is a 1x17 matrix. To fix this, either make A one bigger in each dimension, or make b one shorter. If you make b one shorter, then you need to also make x0 one shorter, because that is 1x17 as well.
n=4;
f=4*ones(1,n^2);
m=diag(f);
format short %Fixed decimal format with a total of 4 digits
j=0;
while j~=n^2
[m(j+1,j+2)]=-1;
[m(j+1,j+4)]=-1;
[m(j+2,j+1)]=-1;
[m(j+4,j+1)]=-1;
j=j+1 ;
end
A=m(1:n^2,1:n^2);
b=zeros(n^2-1, 1);
b(1)=1;
b(n^2)=1;
x0 = (1:n^2)';
L = tril(A,-1);
D = diag(diag(A));
U = -triu(A,1);
K=D+L;
J=K\U;
P=K\b;
x=J*x0+P;
There's no way for me to know if this produces the right output, but it runs without error now.
