y=[0.67 0.67 0.67 0.67 0.67 0.02 0.02 0.02 0.02];
 
temp = 727;
R = 8.3144;
Q = 140000;
temp_k = temp + 273;
d0= 0.000041;
D = d0 * exp( -Q /( R * temp_k)) * 1 * 10^10;
 
dt = 1;
dx = 1;
 
tmp = (D*dt)/dx^2;
 
n = 100;
[s_h s_w] = size(y);
 
y2 = zeros(1, s_w);
 
for j = 1:n
    for i = 1:s_w
 
        im = i-1;
        ip = i+1;
        if im < 1
            im = 1;
        end
        if ip > s_w
            ip = s_w;
        end
        y2(i) = (1 - 2*tmp)*y(i) + tmp*(y(im) + y(ip));
 
    end
 
    y = y2;
end
 
plot(y2, 'DisplayName', 'y2', 'yDataSource', 'y2');

temp = 727;
R = 8.3144;
Q = 140000;
temp_k = temp + 273;
d0= 0.000041;
D = d0 * exp( -Q /( R * temp_k)) * 1 * 10^10;
 
dt = 0.01;
dx = 0.1;
pN = 0.1;
 
if (D*dt)/dx^2 <= 0.5
 
    % Obszar rozwiązania: 100 um    
    %y=[0.67 0.67 0.67 0.67 0.67 0.02 0.02 0.02 0.02];
    y(1:5) = 0.67;
    y(6:100) = 0.02;
 
    ksi = 6;
    for k = 0:dt:10
 
        y = [y(1), y, y(length(y))];
        y1 = y;
 
        for i = 2:1:ksi+1
           y1(i) = D *( dt / dx^2 ) * ( y(i+1) + y(i-1)) + y(i) * (1 - 2*D*(dt/dx^2));
        end
 
        y = y1;
        y = y(2:1:length(y)-1);
 
        plot(y)
        hold on
 
        % rownanie prostej przechodzacej przez dwa punkty
        % (x2 - x1)(y - y1) = (y2 - y1)(x - x1)
        yLine = -0.00407 * temp_k + 3.715;
 
        % zwiekszenie ksi 
        if y(ksi) >= yLine && ksi < length(y)-1
            ksi = ksi + 1;
        end
 
        % zwiekszenie temperatury o zadany prog i kolejna iteracja
        temp_k = temp_k + pN * dt;
        D = d0 * exp( -Q /( R * temp_k)) * 1 * 10^10;
 
    end
   % disp(y);
 
end

temp = 780;
Q = 140000; 
R = 8.3144;
temp_k = temp + 273; 
d0 = 0.000041;
D=d0*exp(-Q /(R * temp_k))*1E10 ;
dt = 0.01; 
dx = 0.1; 
pN = 0;
time = 50;
%Warunek początkowy NEUMANA
if (D*dt)/dx^2<=0.5
 
    K(6:100) = 0.02; %wypełniamy 6-100 elementy stężeniem o wielkości 0.2
    K(1:5) = 0.67; %wypełniamy 1-5 elementy stężeniem o wielkości 0.67
    ksi = 6; %ustawiamy na którym elemencie jest granica ksi
 
    Clen = int32(time/dt); %obliczamy ilość kroków
    C(1:Clen) = 0.0; %tworzymy tablicę wypełnioną zerami do której będzie 
    %obliczane stężenie całego materiału w każdym kroku czasowym
    l = 1; %iterator powyższej tablicy
 
    for k = 0:dt:time %zaczynamy proces
 
        C(l) = trapz(K); %do tablicy zapisujemy stężenie w materiale
        l = l + 1;
 
        K = [K(1), K, K(length(K))]; %zwiększamy rozmiar tablicy z lewej i z prawej strony
        K1 = K; %przepisujemy tablicę do tablicy pomocniczej (k+1)
        K(ksi+1) = K(ksi); %powielamy element graniczny do obliczeń
        for i = 2:1:ksi+1 %obliczamy stężenie do granicy
           K1(i) = D*(dt/dx^2)*(K(i+1)+K(i-1)) + K(i)*(1 - 2*D*(dt/dx^2));
           %obliczanie stężenia dla i-tego elementu
        end
 
        K = K1; %przepisanie tablicy spowrotem
        K = K(2:1:length(K)-1); %zmniejszenie rozmiaru tablicy do stanu początkowego
 
        if (mod(l,100) == 2) %wyświetlamy co 100 element na podstawie iteratora
            plot(K) %wyświetlamy wykres stężenia w materiale
            hold on
        end
 
        y=-0.004074074074074075*temp+3.7155555555555555;
        %na podstawie wyznaczonego wzoru na prostą G-S wyliczamy stężenie
        %dla danej temperatury
        if K(ksi) >= y && ksi < length(K)-1 %jeśli stężenie jest wyższe 
            %niż G-S i ksi nie jest ostatnim elementem to zwiększamy granicę
            ksi = ksi + 1;
        end
 
        temp = temp + pN * dt; %zwiększamy co krok temperaturę
        temp_k = temp + 273;
        D=d0*exp(-Q /(R * temp_k))*1E10; %ponownie obliczamy wsp. dyfuzji
        if (D*dt)/dx^2>0.5 %sprawdzamy czy warunek NEUMANNA nadal jest spełniony
            disp('warunek neumanna przestal byc spelniony dla temp');
            disp(temp);
            break;
        end
    end
    figure;
    plot(C); %wyświetlamy wykres zmiany stężenia co krok
end

temp = 780;
Q = 140000; 
R = 8.3144;
temp_k = temp + 273; 
d0 = 0.000041;
D=d0*exp(-Q /(R * temp_k))*1E10 ;
dt = 0.01; 
dx = 0.1; 
pN = 0;
time = 50;
K(6:100) = 0.02; 
K(1:5) = 0.67; 
ksi = 6; 
 
Clen = int32(time/dt);
C(1:Clen) = 0.0; 
%Neuman
if (D*dt)/dx^2<=0.5
 
    j = 1;
    for k = 0:dt:time
 
        C(l) = trapz(K); %Calka
        j = j + 1;
        K = [K(1), K, K(length(K))];
        TempK = K;
        K(ksi+1) = K(ksi);
        for i = 2:1:ksi+1 
           TempK (i) = D*(dt/dx^2)*(K(i+1)+K(i-1)) + K(i)*(1 - 2*D*(dt/dx^2));
        end
 
        K = TempK;
        K = K(2:1:length(K)-1);
 
        if (mod(j,100) == 2)
            plot(K)
            hold on
        end
 
        y=-0.004074074074074075*temp+3.7155555555555555;
        %zwiekszanie granicy
        if K(ksi) >= y && ksi < length(K)-1 
            ksi = ksi + 1;
        end
 
        temp = temp + pN * dt;
        temp_k = temp + 273;
        D=d0*exp(-Q /(R * temp_k))*1E10;
        if (D*dt)/dx^2>0.5
            break;
        end
    end
    figure;
    plot(C);
end