clear;
clc;
l = 50;
h = 25;
dx = 1;
d = dx;
T_0 = 20;
T_1 = 10;
T_2 = 70;
T_3 = 30;
Kw = 1000;
Kz = 5;
C = 5e4;
Tsym = 20;
dT = 1e-3;
global U;
U = [T_1 T_2 T_3]';
global n m;
n = l/dx;
m = h/dx;
dS = dx*dx;
dC = C/(n*m);
T1_start = round(m/2);
T2_start = round(m/2);
global A B;
A = sparse(n*m, n*m);
B = sparse(n*m, 3);
for i = 1:m
    for j = 1:n
        A_ = sparse(1, n*m);
        B_ = sparse(1,3);
        if i == 1
            B_(3) = B_(3) + Kz*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
        else
            A_(indeks(i-1,j)) = A_(indeks(i-1,j)) + Kw*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
        end
        if i == m
            B_(3) = B_(3) + Kz*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
        else
            A_(indeks(i+1,j)) = A_(indeks(i+1,j)) + Kw*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
        end
        if j == n
            if i < T2_start
                B_(3) = B_(3) + Kz*dS/dC;
                A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
            else
                B_(2) = B_(2) + Kw*dS/dC;
                A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
            end
        else
            A_(indeks(i,j+1)) = A_(indeks(i,j+1)) + Kw*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
        end
        if j == 1
            if i < T1_start
                B_(1) = B_(1) + Kw*dS/dC;
                A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
            else
                B_(3) = B_(3) + Kz*dS/dC;
                A_(indeks(i,j)) = A_(indeks(i,j)) - Kz*dS/dC;
            end
        else
            A_(indeks(i,j-1)) = A_(indeks(i,j-1)) + Kw*dS/dC;
            A_(indeks(i,j)) = A_(indeks(i,j)) - Kw*dS/dC;
        end
        A_(indeks(i,j)) = A_(indeks(i,j)) - 2*Kz*dS/dC;
        B_(3) = B_(3) + 2*Kz*dS/dC;
        A(indeks(i,j),:) = A_;
        B(indeks(i,j),:) = B_;
    end
end
X0 = T_0*ones(n*m,1);
TT = [0 Tsym];
opts = odeset('MaxStep',dT);
[X,T] = ode45('grzanie_plyty',TT,X0);
for i = 1:m
    for j = 1:n
        ZZ(i,j) = T(1,indeks(i,j));
    end
end
h1 = mesh(ZZ);
axis([0 50 0 25 0 70]);
pause;
for k = 1:10:length(T)
    for i = 1:m
        for j = 1:n
            ZZ(i,j) = T(k,indeks(i,j));
        end
    end
    set(h1,'zdata',ZZ);
    pause(0.01);
end

function y=indeks(i,j)
global n;
y=(i-1)*n+j;
end

function Xdot=grzanie_plyty(t,X)
global A B U;
Xdot = A*X + B*U;
end