function r
format long;
m=[1 0 0 0 0;0 1 0 0 0;0 0 1 0 0;0 0 0 1 0;0 0 0 0 1]; 
k=[3 -1 0 0 0;-1 2 -1 0 0;0 -1 2 -1 0;0 0 -1 2 -1;0 0 0 -1 3]; 
n=length(m);
x=zeros(1,n);
z=zeros(n,n);
u=rand(n,1);
d=inv(k)*m;
for i=1:1:n
    [w,y]=diedai(d,m,k,n,u);
    x(1,i)=w;
    z(:,i)=y;
    d=d-(1/w^2)*y*y'*m;
    u=y;
end
w=x
y=z