%Two-dimensional diffusion with the wall between j=m/2 and j=m/2+1
S=[];Sm=[];P=[];q=[];

n=5;m=10;
d=0.2;v=.2;
Sm0=zeros(n,m);Sm0(5,1)=1;
S(:,1)=reshape(Sm0,n*m,1);
%Basic matrix
P=zeros(n*m);
for i=1:n*m-n;P(i,i+n)=d;P(i+n,i)=d+v;end
for i=1:n*m-1;P(i,i+1)=d;P(i+1,i)=d;end
for i=1:m-1;P(i*n,i*n+1)=0;P(i*n+1,i*n)=0;end;
%creating the wall
for i=1:n;P(n*m/2+i,n*m/2-n+i)=0;end;P(n*m/2+1,n*m/2-n+1)=d+v;
for i=1:n;P(n*m/2-n+i,n*m/2+i)=0;end;P(n*m/2-n+1,n*m/2+1)=d;
%Diagonal elements
for i=1:n*m;P(i,i)=1-sum(P(:,i));end;

N=260;%Number of transitions
for k=1:N
   S(:,k+1)=P*S(:,k);
end
figure(1);colormap(autumn)
%VARIANTS OF GRAPH PRESENTATION. TRY ALL.

for i=1:N;Sm=reshape(S(:,i),n,m);bar3(Sm');axis([.5 n+.5 .5 m+.5 0 .25*max(max(Sm0))]);pause(.1);hold off;end

%for i=1:N;Sm=reshape(S(:,i),n,m);surf(Sm);axis([.5 m+.5 .5 n+.5 0 .25*max(max(Sm0))]);pause(.1);hold off;end

%for i=1:N;Sm=reshape(S(:,i),n,m);contourf(Sm,20);axis([.5 m+.5 .5 n+.5]);pause(.1);hold off;end

%Sm=reshape(S(:,N),n,m);bar3(Sm);axis([.5 m+.5 .5 n+.5 0 max(max(Sm0))]);hold off
