% Smisena uloha pro vlnovou rovnici         
% du/dt = p d2u/dx2 + f(x,t) ... v (a,b)x(0,T)             
%                       
% pocatecni: u(x,0) = phi(x) ... v <a,b>              
% okrajove:          
% u(a,t) = alfa(t) v <0,T>              
% u(b,t) = beta(t) v <0,T>    
%
%vstupy
p = 3.2;               % difuzni koeficient, kladna konstanta,
f = @(x,t)(4*x + 2*t);     % funkce vnejsich zdroju,

a = 3;   % interval promenne x
b = 11;

phi   = @(x)-1+11*x-x.^2;   % pocatecni podminka phi(x)

G=phi(a)
H=phi(b)

alpha = @(t) G + log(1+t);      % okrajova podminka alpha(t)
beta  = @(t) H + log(1+t);      % okrajova podminka beta(t)

h   = 2;   % prostorovy krok      
tau = 0.5;   % casovy krok

% Overeni podminek souhlasu
 if phi(a)==alpha(0) && phi(b)== beta(0)
    disp('Podminky souhlasu jsou splneny.');
 else
    disp('Podminky souhlasu nejsou splneny.');
    return;
 end   
sigma = p * tau / h^2;

% Overeni podminky stability
if sigma <= 0.5 
   v=['Podminka stability pro h=',num2str(h),...
      ', tau=',num2str(tau),...
      ' a  p=',num2str(p),' je splnena. sigma = ',num2str(sigma)];
 disp(v);
else
   v=['Podminka stability pro h=',num2str(h),...
         ', tau=',num2str(tau),...
      ' a p =',num2str(p),' neni splnena! sigma = ',num2str(sigma)];
 disp(v);
end

% Vypocet 
% body xi
x = a : h :b+h/10 %ověřit, zda x(end) je opravdu b
pocet_bodu=length(x);
% nulta vrstva
U0 = phi(x)
% EXPLICITNI SCHEMA
% prvni vrstva
U1E(1) = alpha(tau);
U1E(pocet_bodu) = beta(tau);
for i = 2:pocet_bodu-1
    U1E(i)=sigma*(U0(i-1)+U0(i+1))+(1-2*sigma)*U0(i)+tau*f(x(i),0);
end
U1E
% druha vrstva
U2E(1) = alpha(2*tau);
U2E(pocet_bodu) = beta(2*tau);
for i = 2:pocet_bodu-1
    U2E(i)=sigma*(U1E(i-1)+U1E(i+1))+(1-2*sigma)*U1E(i)+tau*f(x(i),tau);
end
U2E
% % treti vrstva
% U3E(1) = alpha(3*tau);
% U3E(pocet_bodu) = beta(3*tau);
% for i = 2:pocet_bodu-1
%     U3E(i)=sigma*(U2E(i-1)+U2E(i+1))+(1-2*sigma)*U2E(i)+tau*f(x(i),2*tau);
% end

% IMPLICITNI SCHEMA: soustava rovnic pro prvni casovou vrstvu
% rovnice pro vnitrni body,rozmer matice (pocet_vnitrnich x pocet_vnitrnich)
pocet_vnitrnich = pocet_bodu-2;
% matice soustavy
Uimplicitni = zeros(pocet_vnitrnich,pocet_vnitrnich); %i když je 3-diagonalni
for i=1:pocet_vnitrnich Uimplicitni(i,i) = (1+2*sigma); end %hlavni diagonala
for i=2:pocet_vnitrnich Uimplicitni(i,i-1) = -sigma; end    %dolni diagonala
for i=1:pocet_vnitrnich-1 Uimplicitni(i,i+1) = -sigma; end  %horni diagonala

% prava strana
prava_strana = tau * f(x(2:end-1),tau) + U0(2:end-1);
prava_strana(1) = prava_strana(1) + sigma*alpha(tau);
prava_strana(end) = prava_strana(end) + sigma*beta(tau);

%
Uimplicitni
prava_strana



