% Grey-Scott egyenlet 2D-ben, stabil mintáztaképzõdés
    %du/dt = Du*du^2/dx^2-u*v^2+F*(1-u);
    %dv/dt = Dv*du^2/dx^2+u*v^2-(C+F)*v;

dt = 0.2; % idõbeli lépésköz %D*2*dt/dx^2<1
dx = 1/100; %térbeli lépésköz
dy = 1/100; %térbeli lépésköz

x = 0:dx:1;
y = 0:dy:1;
nx = length(x);
ny = length(y);
nt = 20000; %lépésszám

% diffúziós együtthatók: diff sebességét határozza meg
Du=6e-5;
Dv=2e-5;

% paraméterek:
%F=0.06 sávok, 0.03 foltok végül
C = 0.065;%0.015; %fogyás 
F = 0.06;%0.03; %termelõdés 

% kezdeti feltételek: u=1; v=0 alapból (0-1 közt cc-k) 
%más kezdeti feltételeket is nézni, mikortól lesz mintázat?
u = ones(nx,ny);
%u(20:30, 10:40) = 0; 
v = zeros(nx,ny);
v(20:40, 20:40) = 1; %kezdeti állapot (négyzet)

%ábrázolás: u és v mátrixok ábrázolása
figure;
subplot(1,2,1)
    imagesc(u);
     title(['lépés: 0   u']);
    caxis([0 1]); %lecsengés normált színekkel jobban látszódik
subplot(1,2,2);
    imagesc(v);
    title(['v']);
    caxis([0 1]);
    
   waitforbuttonpress;
   
for t = 1:nt %idõlépések
    
    %
    %mátrix mûveletekkel: sokkal "matlabosabb" és gyorsabb
    u_iprev=[u(end, :); u(1:end-1, :)]; % (u2(i-1,j) helyett) mátrix, ami minden pozícióban az elõzõ sorban lévõ elemet tartalmazza +periodikus határ 
    u_iaft=[u(2:end, :); u(1,:)];
    u_jprev=[u(:, end), u(:, 1:end-1)];
    u_jaft=[u(:, 2:end), u(:, 1)];
    
    v_iprev=[v(end, :); v(1:end-1, :)]; % (u2(i-1,j) helyett) mátrix, ami minden pozícióban az elõzõ sorban lévõ elemet tartalmazza +periodikus határ 
    v_iaft=[v(2:end, :); v(1,:)];
    v_jprev=[v(:, end), v(:, 1:end-1)];
    v_jaft=[v(:, 2:end), v(:, 1)];
    
    u= u+dt*(-u.*v.^2+F*(1-u)+Du/(dx^2)*(u_iprev+u_jprev-4*u+u_iaft+u_jaft));
    v= v+dt*(+u.*v.^2-(C+F)*v+Dv/(dx^2)*(v_iprev+v_jprev-4*v+v_iaft+v_jaft));
    %}
    
    %{
    u2=u;
    v2=v;
    for i =1:nx   
        for j =1:ny 
            
            if i==1 %periodikus peremek
               i_before=nx;
            else 
                i_before=i-1;
            end
            if i==nx
                i_after=1;
            else
                i_after=i+1;
            end
            
            if j==1
               j_before=ny;
            else 
                j_before=j-1;
            end
            if j==ny
                j_after=1;
            else
                j_after=j+1;
            end
            
             u(i, j) = u2(i, j) + dt*(-u2(i,j))*v2(i,j)^2 +dt*F*(1-u2(i,j)) +...
                                 Du*dt/dx^2*1*(u2(i_after, j)-2*u2(i, j)+u2(i_before, j))...
                                + Du*dt/dy^2*1*(u2(i, j_after)-2*u2(i, j)+u2(i, j_before));
            
             v(i, j) = v2(i, j) + dt*u2(i,j)*v2(i,j)^2 -(C+F)*dt*v2(i,j) +...
                                 Dv*dt/dx^2*1*(v2(i_after, j)-2*v2(i, j)+v2(i_before, j))...
                                + Dv*dt/dy^2*1*(v2(i, j_after)-2*v2(i, j)+v2(i, j_before));                            
  
        end     
    end
    %}
    
    if mod(t,50)==0 %minden x. lépésben rajzolok csak ki
        
    subplot(1,2,1)
    imagesc(u); 
    caxis([0 1]); %lecsengés normált színekkel jobban látszódik
    title(['lépés: ', num2str(t), ' u']);
    subplot(1,2,2);
    imagesc(v);
    title(['v']);
    caxis([0 1]);
    
   %waitforbuttonpress;
    pause(0.001);
    end
    
end