clear;
close all;
%------------------------------
% Solution of the 2d RD equation by use of ADI-method
%---------------------------------
%-----------------------------
du=5;
dv=12;
lam=9;
k=3;
homu=k;
homv=lam/k;
%-----------------------------
L=50; % length 
M=100; % amount of points on the (0,L)
x=linspace(0, L, M);
dx=x(2)-x(1);% space grid size
dt=0.05; % time grid size
a=(du*dt)/(2*dx^2);
b=(dv*dt)/(2*dx^2);
T=30000; % amount of time steps
tt=2000;
rad=0.01;
nn=1;
%-----------------------------
a1=-a;
a2=1+2*a;
a3=-a;
%-----------------------------
b1=-b;
b2=1+2*b;
b3=-b;
%-----------------------------
[x,y]=meshgrid(0:dx:L);

    rausch=0.1*(rand(M,M)-0.5);
    u0=homu+rausch;
    v0=homv+rausch;

    u=u0;
    v=v0;
    
[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')
saveas(gcf,['bild_' num2str(0) '.png'])
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;	

triB = spdiags(ones(M,1)*[-b 1+2*b -b],-1:1,M,M);
triB(1,2)=-2*b;						
triB(M,M-1)=-2*b;

A = spdiags(ones(M,1)*[a 1-2*a a],-1:1,M,M);
A(1,2)=2*a;						
A(M,M-1)=2*a;

B = spdiags(ones(M,1)*[b 1-2*b b],-1:1,M,M);
B(1,2)=2*b;						
B(M,M-1)=2*b;

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

 for t=1:T
  % first half step: implicit in x-direction
  u1=u;
  v1=v;
  for i=1:M
  rhsu=A*u1(i,:)'+(dt/2)*((u1(i,:).^2)'.*v1(i,:)'-(lam+1)*u1(i,:)'+k);
  rhsv=B*v1(i,:)'+(dt/2)*(-(u1(i,:).^2)'.*v1(i,:)'+lam*u1(i,:)');
  u(i,:)=rhsu'/triA';
  v(i,:)=rhsv'/triB';
  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
  u1=u;
  v1=v;
  for i=1:M
      rhsu=A*u1(:,i);%+(dt/2)*(lam*u1(:,i)-(u1(:,i).^3)-v1(:,i)+k);
      rhsv=B*v1(:,i);%+(dt/2)*(u1(:,i)-v1(:,i));
      u(:,i)=triA\rhsu;
      v(:,i)=triB\rhsv;
  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')
       saveas(gcf,['bild_' num2str(t) '.png'])
       hold on 
       pause(0.5)
       hold off
   end;
 end;
 