clear;
close all;
%------------------------------
% Solution of the 2d heat equation by use of ADI-method
%---------------------------------

L=1; % length of the string
M=100; % amount of points on the (0,L)
D=1;
T=20;
x=linspace(0, L, M);
y=linspace(0, L, M);
dx=x(2)-x(1);% space grid size
dt=0.001;
a=(D*dt)/(2*dx^2);
ys=1;
tt=5;
%-----------------------------
[x,y]=meshgrid(0:dx:L);
f=inline('exp(-(20*(z-0.5)).^2-(20*(k-0.5)).^2)*1','z','k'); % an initial wave u(x,)
u=f(x,y);

[C,h]=contourf(x,y,u,50);
axis square
set(h,'LineStyle','none')
colorbar
colormap bone
title(['t=' num2str(0)],'Color','b')
xlabel('x')
ylabel('y')
clear h
hold on 
pause(0.5)
hold off
% % %-----------------------------
triA = spdiags(ones(M,1)*[-a 1+2*a -a],-1:1,M,M);
triA(1,2)=-2*a;						
triA(M,M-1)=-2*a;	
%matrix A
%  maina=(1+2*a)*ones(M,1);
%  an=-a*ones(M-2,1);
%  an(M-1)=-2*a;
%  av=-a*ones(M-1,1);
%  av(1)=-2*a;
%  triA=diag(maina)+diag(an,-1)+diag(av,1);
% %matrix B
 mainb=(1-2*a)*ones(M,1);
 bn=a*ones(M-2,1);
 bn(M-1)=2*a;
 bv=a*ones(M-1,1);
 bv(1)=2*a;
 triB=diag(mainb)+diag(bn,-1)+diag(bv,1);

 % %-----------------------------

 for t=1:T
  % first half step: implicit in x-direction
  v=u;
  for i=1:M
  rhs=triB*v(i,:)';
  u(i,:)=rhs'/triA';
  end;
%   %i=1;
%   rhs=(1-2*a)*v(1,:)+2*a*v(2,:);
%   u(1,:)=rhs/triA';
%   for i=2:M-1
%       rhs=a*v(i-1,:)+(1-2*a)*v(i,:)+a*v(i+1,:);
%       u(i,:)=rhs/triA';
%   end;
%   %i=M;
%   rhs=(1-2*a)*v(M,:)+2*a*v(M-1,:);
%   u(M,:)=rhs/triA';
  %-----------------------------  
  % second half-step: implicit in y-direction
  v=u;
  for i=1:M
      rhs=triB*v(:,i);
      u(:,i)=triA\rhs;
  end;
%   %i=1;
%   rhs=(1-2*a)*v(:,1) +2*a*v(:,2);	
%   u(:,1) = triA\rhs;
%   for i=2:M-1
%       rhs=a*v(:,i-1)+(1-2*a)*v(:,i)+a*v(:,i+1);
%       u(:,i)=triA\rhs;
%   end;
%   %i=M;
%   rhs=(1-2*a)*v(:,M)+2*a*v(:,M-1);
%   u(:,M)=triA\rhs;
   %----------------------------- 
    if mod(t,tt)==0
      [C,h]=contourf(x,y,u,50);
       axis square
       set(h,'LineStyle','none')
       colorbar
       colormap gray
       title(['t=' num2str(t)],'Color','b')
       xlabel('x')
       ylabel('y')
       hold on 
       pause(0.5)
       hold off
   end;
 end;
 
 

%  for t=1:T % time loop
%   % first half step: implicit in x-direction
%   v=u; % old values
%   for i=1:M
%       rhs=triB*v(i,:)';
%       u(i,:)=inv(triA)*rhs;
%   end; 
%   v=u;
%   % second half-step: implicit in y-direction
%   for i=1:M
%       rhs1=triB*u(:,i);
%       u(:,i)=inv(triA)*rhs1;
%   end;
%   if mod(t,tt)==0
%       [C,h]=contourf(x,y,u,50);
%        axis square
%        set(h,'LineStyle','none')
%        colorbar
%        title(['t=' num2str(t)],'Color','b')
%        xlabel('x')
%        ylabel('y')
%        hold on 
%        pause(0.5)
%        hold off
%        clear h
%    end;
%  end;
%   