function [U,t] = THETA(eq,t,U0,th,log,num_impl)
%THETA Řeší nestacionární problém vedení tepla theta metodu
% Pro theta=0.5 se jedná o Crankovu-Nicolsonovu metodu. Při nesouladu
% počátečních a okrajových podmínek je vhodné nastavit parametr num_impl
% na malé číslo např. 5. Pak se prvních 5 kroků bude počítat implicitní
% Eulerovou metodou, která sice není tak přesná, je ovšem stabilní.
% Vstupy: eq - strukture představující data přímé úlohy vedení tepla, t -
% dělení časového intervalu, U0 - počáteční podmínka v uzlech triangulace,
% th - theta viz theta metoda, log - má-li být použita logaritmická časová
% škála, num_impl - počet kroků, kdy má být použita implicitní Eulerova
% metoda
% Výstupy: U - matice řešení, t - dělení časového intervalu
t = t(:);
t = t';
if eq.stationary
    error('Úloha je stacionární!');
end
% Časové kroky jsou buď dodány, nebo si je vygeneruji
if isempty(t)
    % Zejména při nesouladu okrajových a počátečních podmínek je vhodné
    % použít logaritmické dělení časového interval
    if log
        t = logspace(-3,log10(eq.time(2)),eq.num_tsteps-1);
        t = [eq.time(1),t];
        if t(end) ~= eq.time(2)
            t(end) = eq.time(2);
        end
        
    else
        t = linspace(eq.time(1),eq.time(2),eq.num_tsteps);
    end
end
% Inicializace matice, kam bude uloženo řešení
U = zeros(size(eq.phi(eq.domain.p(1,:)',eq.domain.p(2,:)),1)',size(t,2));
% Počáteční podmínka je zadaná buď jako funkce, nebo v každém uzlu triangulace
if isempty(U0)
    U(:,1) = eq.phi(eq.domain.p(1,:),eq.domain.p(2,:));
else
    U(:,1) = U0;
end
% Kontrola parametru theta
if th < 0.5
    warning('Pozor, theta metoda s th < 0.5 může být pro velké časové kroky nestabilní!');
end
if th > 1 || th < 0
    error('Parametr th musí být z intervalu <0,1>!');
end
warning('off');
% Hodnotu parametru theta si dočasně uložím
th_last = th;
% Nastavené fsolve parametrů
options = optimoptions('fsolve','Display','iter','Algorithm','levenberg-marquardt');
% Theta metoda
if eq.linear
    % Lineární úloha
    for i = 1:size(t,2)-1
        % Implicitní Eulerova metoda po num_impl kroků
        if i < num_impl
            th = 1;
        else
            th = th_last;
        end
        tau_i = t(i+1)-t(i);
        t_ith = t(i)+th*tau_i;
        [C,K,F] = FEM(eq,t_ith);
        U(:,i+1) = (C+tau_i*th*K)\((C-tau_i*(1-th)*K)*U(:,i)+tau_i*F);
    end
else
    % Nelineární úloha
    for i = 1:size(t,2)-1
        % Implicitní Eulerova metoda po num_impl kroků
        if i < num_impl
            th = 1;
        else
            th = th_last;
        end
        tau_i = t(i+1)-t(i);
        t_ith = t(i)+th*tau_i;
        [U(:,i+1),~,exitflag] = fsolve(@(U_act)wrapper(U_act,eq,tau_i,th,t_ith,U(:,i)),U(:,i),options);
        if exitflag < 0
            error('Nelineární soustava rovnic nebyla vyřešena. Zastavuji theta metodu.');
        end
    end
end
warning('on');
end

function R = wrapper(U_act,eq,tau_i,th,t_ith,U_prev)
    U_t_ith = U_prev+th*(U_act-U_prev);
    [C,K,F] = FEM(eq,t_ith,U_t_ith);
    R = (C+tau_i*th*K)*U_act-((C-tau_i*(1-th)*K)*U_prev+tau_i*F);
end

