clear;
clc;
Kw = [1000 2000 1000];
Kz = [1 10 1];
Kz1 = 1; Kz2 = 10; Kz3 = 1;
Kw1 = 1000; Kw2 = 2000; Kw3 = 1000;
T0 = 20;
T1 = 10;
T2 = 70;
T3 = 30;
global A B U;
U = [T1 T2 T3]';
C = 5e4;
l = 15;
dx = 0.1;
dS = dx^2;
n = l/dx;
dC = C/n;
A = sparse(n,n);
B = sparse(n,3);
Xstart = T0*ones(n,1);
for i=1:n
    A_ = sparse(1,n);
    B_ = sparse(1,3);
    if i<=50
        A_(i) = A_(i) - 4*Kz1*dS/dC;
        B_(3) = B_(3) + 4*Kz1*dS/dC;
        if i ==1
            A_(i) = A_(i) - (Kw1*dS/dC);
            B_(1) = B_(1) + (Kw1*dS/dC);
        else
            A_(i) = A_(i) - (Kw1*dS/dC);
            A_(i-1) = A_(i-1) + (Kw1*dS/dC);
        end
        if i == 50
            A_(i) = A_(i) - (((Kw1+Kw2)/2)*dS/dC);
            A_(i+1) = A_(i+1) + (((Kw1+Kw2)/2)*dS/dC);
        else
            A_(i) = A_(i) - (Kw1*dS/dC);
            A_(i+1) = A_(i+1) + (Kw1*dS/dC);
        end
    elseif i>50
        if i<=100
            A_(i) = A_(i) - (4*Kz2*dS/dC);
            B_(3) = B_(3) + (4*Kz2*dS/dC);
            if i == 50
                A_(i) = A_(i) - (((Kw1+Kw2)/2)*dS/dC);
                A_(i-1) = A_(i-1) + (((Kw1+Kw2)/2)*dS/dC);
            else
                A_(i) = A_(i) - (Kw2*dS/dC);
                A_(i-1) = A_(i-1) + (Kw2*dS/dC);
            end
            if i == 100
                A_(i) = A_(i) - (((Kw3+Kw2)/2)*dS/dC);
                A_(i+1) = A_(i+1) + (((Kw3+Kw2)/2)*dS/dC);
            else
                A_(i) = A_(i) - (Kw2*dS/dC);
                A_(i+1) = A_(i+1) + (Kw2*dS/dC);
            end
        else
            A_(i) = A_(i) - (4*Kz3*dS/dC);
            B_(3) = B_(3) + (4*Kz3*dS/dC);
            if i == 100
                A_(i) = A_(i) - (((Kw3+Kw2)/2)*dS/dC);
                A_(i-1) = A_(i-1) + (((Kw3+Kw2)/2)*dS/dC);
            else
                A_(i) = A_(i) - (Kw3*dS/dC);
                A_(i-1) = A_(i-1) + (Kw3*dS/dC);
            end
            if i == 150
                A_(i) = A_(i) - (Kw3*dS/dC);
                B_(2) = B_(2) + (Kw3*dS/dC);
            else
                A_(i) = A_(i) - (Kw3*dS/dC);
                A_(i+1) = A_(i+1) + (Kw3*dS/dC);
            end
        end
    end
    A(i,:) = A_;
    B(i,:) = B_;
end
Tsym = 50000;
[t,X] = ode45('grzanie_preta2', [0 Tsym], Xstart);
h = line(1:length(X(1,:)), X(1,:),'erasemode', 'xor');
axis([0 n 0 max(U)]);
pause;
for i=2:10:length(t);
    set(h, 'ydata', X(i,:));
    drawnow;
    pause(0.015);
end;

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