% 1. gyakorlat: diffúzió Dirichlet és Neumann peremmel

dt = 0.1; %idõbeli lépésköz
dx = 1; %térbeli lépésköz
D=4;%1;%4;%5;%6;% %diffúzió sebessége %ezt változtassuk most
nt = 500; %idõtartomány vége, eddig futtatok

x = 1:dx:200; %viszgált pontok egy egyenesen
nx = length(x); %vizsgált pontok száma
u = zeros(nx,1); % kezdeti értékek 0-ák
U=zeros(nx, nt); %ha az idõbeli változásokat 1 ábrán akarjuk nézni

%u(95:104)=1; %csúcs indul, sûrûség függvény,
%u(170:180)=1;
%u(140:150)=1; %két csúcs induláskor
u=0.4+sin(1/5*x)/4; % más kezdeti érték 

figure;

for t = 1:nt %minden idõpillanatban
    u2 = u; %u2=értékek az elõzõ idõplillanatban; u=új értékek
    for i = 1:nx %minden pontban kiszámolom a hullám aktuális értékét   
       %minden pontra kiszámoljuk az új értékét:               
       if i==1 %peremek
             u(i)=0; %1. fix perem (Dirichlet perem): alapból is ez volt, torzít a széleken, lefolyik az anyag, 0 lesz
             % u(i) = u2(i) + D*dt/dx^2*(u2(i+1)-2*u2(i)+u2(end)); % 3. periodikus határ: az elsõ és az utolsó pozíció összekapcsolva, kör lesz, kiegyenlítõdik az átlaghoz 
              u(i) = u2(i) + D*dt/dx^2*(u2(i+1)-1*u2(i)); %2. Neumann perem: gradiens 0 széleken=utolsó és utolsó elötti mindig ugyanazon az értéken van, kiegyenlítõdik az átlaghoz
       elseif i>1 && i<nx
             u(i) = u2(i) + D*dt/dx^2*(u2(i+1)-2*u2(i)+u2(i-1)); %dU/dt=D*dU^2/dx^2 %diffúzió, szétfolyik
       else %peremek
            u(i)=0; %1. fix perem (Dirichlet perem): alapból is ez volt, torzít a széleken felolyik az anyag
           u(i) = u2(i) + D*dt/dx^2*(-1*u2(i)+u2(i-1)); %2. Neumann perem: gradiens 0 széleken=utolsó és utolsó elötti mindig ugyanazon az értéken van, kiegyenlítõdik az átlaghoz
          % u(i) = u2(i) + D*dt/dx^2*(u2(1)-2*u2(i)+u2(i-1)); % 3. periodikus határ: az elsõ és az utolsó pozíció összekapcsolva, kör lesz, kiegyenlítõdik az átlaghoz

       end
    end
    
    %if mod(t,100)==0 %hosszú futásnál minden 100.-at rajzol csak ki
    plot(u); %aktuális értékeket kirajzol
    grid(gca, 'minor'); % sûrû rácsot tesz a képre
    ylim([-0.5 1.5]); %ábrázolás határai
    xlim([1, 200]);%xlim([10, 190]);%
    title(['lépés: ', num2str(t)]); %kiírítja az ábrára, hogy hányadik lépésnél vagyok
    %waitforbuttonpress; %gombnyomásra lép egyet a ciklus
     pause(0.01); %ilyen idõközönként rajzolja ki az új állapotot
     %end
     
    U(:, t)=u;
end

mesh(U); %ha az idõbeli változásokat 1 ábrán akarjuk nézni

%halad a csúcs ha u(i+1)=u2(i) + D*dt/dx^2*(u2(i+1)-2*u2(i)+u2(i-1));
% a rendszer elszáll, ha:
%D*2*dt/dx^2>1