function [u,p,k]=SplitIsotropic2(f,X,Tol, lambda,mu)
tic
n=size(f,1); m=size(f,2);
n1=n-1; m1=m-1;
u=f; dx =zeros(n,m); dy=zeros(n,m); bx=zeros(n,m); by=zeros(n,m);
k=0; p(1) = norm(f-X,2)/norm(X,2);nn(1) = norm(u,2);
while nn(k+1) >Tol
k=k+1;
u1 = u ;
for i=2:n1
for j=2:m1
u(i,j)= g(u,dx,dy, bx, by,f, lambda, mu, i,j);% zhong jian
end
end
%Compute the sx¡¯s
for i=2:n1
for j=2:m1
s(i,j) = sqrt( abs( (u(i+1,j)-u(i-1,j))/2 + bx(i,j))^2 ...
+ abs((u(i,j+1)-u(i,j-1))/2 + by(i,j))^2);%zhong jian
end
end
for i=1:n1
j=1;
s(i,j) = sqrt( abs( (u(i+1,j)-u(i,j)) + bx(i,j))^2 ...
+ abs((u(i,j+1)-u(i,j)) + by(i,j))^2);%xia
end
for j=1:m1
i=1;
s(i,j) = sqrt( abs( (u(i+1,j)-u(i,j)) + bx(i,j))^2 ...
+ abs((u(i,j+1)-u(i,j)) + by(i,j))^2);%zuo
end
for i=2:n
j=m;
s(i,j) = sqrt( abs( (u(i,j)-u(i-1,j)) + bx(i,j))^2 ...
+ abs((u(i,j)-u(i,j-1)) + by(i,j))^2);%shang
end
for j=2:m
i=n;
s(i,j) = sqrt( abs( (u(i,j)-u(i-1,j)) + bx(i,j))^2 ...
+ abs((u(i,j)-u(i,j-1)) + by(i,j))^2);%you
end
for j=1:m1
i=1;
s(i,j) = sqrt( abs( (u(i+1,j)-u(i,j)) + bx(i,j))^2 ...
+ abs((u(i,j+1)-u(i,j)) + by(i,j))^2);%zuo
end
%Compute the d¡¯s
for i=2:n1
for j=1:m
dx(i,j)= (s(i,j)*lambda*((u(i+1,j)-u(i-1,j))/2+bx(i,j)))...
/(s(i,j)*lambda + 1);%zhong jian
end
end
for j=1:m
i=1;
dx(i,j)= (s(i,j)*lambda*((u(i+1,j)-u(i,j))+bx(i,j)))...
/(s(i,j)*lambda + 1);%zuo
end
for j=1:m
i=n;
dx(i,j)= (s(i,j)*lambda*((u(i,j)-u(i-1,j))+bx(i,j)))...
/(s(i,j)*lambda + 1);%you
end
for i=1:n
for j=2:m1
dy(i,j) = (s(i,j)*lambda*(((u(i,j+1)-u(i,j-1))/2 + by(i,j))))...
/(s(i,j)*lambda +1 );%zhong jian
end
end
for i=1:n
j=1;
dy(i,j) = (s(i,j)*lambda*(((u(i,j+1)-u(i,j)) + by(i,j))))...
/(s(i,j)*lambda +1 );%xia
end
for i=1:n
j=m;
dy(i,j) = (s(i,j)*lambda*(((u(i,j)-u(i,j-1))/2 + by(i,j))))...
/(s(i,j)*lambda +1 );%shang
end
%Compute the b¡¯s
for i=2:n1
for j=1:m
bx(i,j)= bx(i,j) + ((u(i+1,j)-u(i-1,j))/2 - dx(i,j));%zhong jian
end
end
for j=1:m
i=1;
bx(i,j)= bx(i,j) + ((u(i+1,j)-u(i,j)) - dx(i,j));%zuo
end
for j=1:m
i=n;
bx(i,j)= bx(i,j) + ((u(i,j)-u(i-1,j)) - dx(i,j));%you
end
for i=1:n
for j=2:m1
by(i,j) = by(i,j) + ((u(i,j+1)-u(i,j-1))/2 - dy(i,j));%zhong jian
end
end
for i=1:n
j=1;
by(i,j) = by(i,j) + ((u(i,j+1)-u(i,j)) - dy(i,j));%xia
end
for i=1:n
j=m;
by(i,j) = by(i,j) + ((u(i,j)-u(i,j-1)) - dy(i,j));%shang
end
nn(k+1) = norm(u-u1,2)/norm(u,2);
p(k+1)= norm(u-X,2)/norm(X,2);
end
toc
end
