function prikladNestacionarni_1d()

% c*du/dt - d/dx (k du/dx(x)) + u = 10+1 v bodech x z (0,1) a casech t z (0,T)
% -k(0,t) du(0,t)/dx = 1*(1-u(0,t)) pro vsechny casy t z <0,T>
% -k(1,t) du(1,t)/dx = 1*(u(1,t)-1) pro vsechny casy t z <0,T>
% u(x,0) = u0(x) v bodech x z [0,1]

l = 1; q = 1; alfa = 1; f = 10; ustriska = 1;
kvlevo = 1; kvpravo = 4;
c = 1;
u0 = ustriska;

ft = @(t) (cos(t));

m = 32;
n = m+1; h = 1/m;

% cas. diskretizace
T = 10; N = 2000;
t = linspace(0,T,N);
dt = T/(N-1);

f = (f+ustriska*q)*ones(1,m);
k = [kvlevo*ones(1,m/2),kvpravo*ones(1,m/2)];
q = q*ones(1,m);
cdt = c/dt * ones(1,m);

% diskretizace
x = linspace(0,l,n);
E = [1:(n-1);2:n];

% okr. podm.
Dirichlet = []; uD = [];
Neumann = [1,n]; tau = alfa*ustriska*[1,1]';

% MKP
Ae1 = @(ke,xe) ( ke/(xe(2)-xe(1))*[1,-1;-1,1] );
Ae2 = @(qe,xe) ( qe*(xe(2)-xe(1))*[1/3,1/6;1/6,1/3] );
be = @(fe,xe) (fe*(xe(2)-xe(1))/2*[1;1]);
[A,b] = MKP(x,E,{Ae1,Ae2},{k,q},{be},{f},Dirichlet,uD,Neumann,tau);
% modifikace A
A(1,1) = A(1,1) + alfa;
A(n,n) = A(n,n) + alfa;
% matice pro casovou diskretizaci
Me = @(cdte,xe) ( cdte*(xe(2)-xe(1))*[1/3,1/6;1/6,1/3] );
M = MKP(x,E,{Me},{cdt},{},{},Dirichlet,uD,[],[]);

% stacionarni reseni
u = A\b;
plot(x,u','r')

uk = u0*ones(n,1);
for k=1:N
    tk = k*dt;
%     % dopredny Euler
%     rhsk = b - A*uk + M*uk;
%     uk = M \ rhsk;
    % zpetny Euler
    rhsk = b*ft(tk) + M*uk;
    uk = (M+A) \ rhsk;
    plot(x,u','r',x,uk','b')
    pause(0.1)
end
