% This demo shows behaviour of stiff ODEs
% We solve the test equation
%   y'(t) = A y(t)
%
% for A random diagonizable matrix with
% eigenvalues
%    lambda_1 = a_1 + i b_1
%    lambda_2 = a_2 + i b_2
%    lambda_3 = lambda_1*
%    lambda_4 = lambda_2*
%
% more detailed:
% let Q be a random orthogonal real matrix
% A = Q'*C*Q,
% with C= [  a_1 b_1
%           -b_1 a_1
%                     a_2 b_2
%                    -b_2 a_2 ]

clear all;

% Set problem data
T      = 5;
tspan  = [0,T];

% eigenvalue parameters
a_1    = -1000;
b_1    = 1000;
a_2    = -1;
b_2    = 6;

a_1    = -2;
b_1    = 0;
a_2    = -2;
b_2    = 0;

% C= matrix with these eigenvalues
C= [a_1,b_1,0,0;...
    -b_1,a_1,0,0;...
    0,0,a_2,b_2;...
    0,0,-b_2,a_2],


% generate "random" orthogonal matrix
%Qraw=randn(4,4);
Qraw =[
    -0.7258    1.3419   -1.4517    0.1602
    -0.8665   -0.9884   -0.6187    0.2874
    -0.4218    1.8179    0.9345    0.6329
    -0.9427   -0.3744    1.0559   -1.4590];
[Q,bla]=qr(Qraw);

% A= "random" matrix with the specified eigenvalues
A=Q'*C*Q;

%eig_C=eig(C),
%eig_A=eig(A),

% "random" initial values
%y0     = randn(4,1);
y0 =[-0.3479;1.2901;1.3412;-0.5808];

% save figures?  (0 not, 1 save)
savefigs=0;

% Set method
%for themethod=1:4
for themethod=2
    
    switch themethod
        case 1
            exp_imp='exp';
            grid='uniform';
        case 2
            exp_imp='imp';
            grid='uniform';
        case 3
            exp_imp='exp';
            grid='switched';
        case 4
            exp_imp='imp';
            grid='switched';
    end
    
    
    % define exact solution
    sol = @(t) expm(A*t)*y0;
    
    % Set list of trial step numbers
    nlist = logspace(1,4,10)%5,40);
    
    % Prepare the figures
    figure(1); clf; hold off
    figure(2); clf; hold off
    
    count=0;
    
    % Loop over all trial step numbers
    for N = nlist
        
        % define the grid
        switch grid
            case 'uniform'
                t = linspace(0,T,2*N);
            case 'switched'
                ts=-log(N)/(a_1);
                t1=linspace(0,ts,N);
                t2=linspace(ts,T,N);
                t=[t1,t2(2:end)]';
        end
        
        nt = length(t);
        ny = length(y0);
        y = zeros(nt,ny);
        y(1,:) = y0';
        
        % calculate the solution
        switch exp_imp
            case 'imp'
                % implicit Euler
                for i=1:nt-1
                    hloc=t(i+1)-t(i);
                    y(i+1,:) = ( (eye(ny)-hloc*A) \ y(i,:)' )';
                end
            case 'exp'
                % explicit Euler
                for i=1:nt-1
                    hloc=t(i+1)-t(i);
                    y(i+1,:) = y(i,:) + (hloc*A *y(i,:)' )';
                end
            otherwise
                exp_imp
                error('exp_imp switch failed');
        end
        % only for error plot define an equivalent h
        h=1/(2*N);
        
        
        % Compute the error (discrete maximum norm on the grid)
        sol_t=zeros(nt,ny);
        E=0;
        for ti=1:nt;
            sol_t(ti,:)= sol(t(ti));
            %E = max(abs(y - sol(t)));
            E=max([E,norm(y(ti,:)-sol_t(ti,:))]);
        end
        sol_t
        
        
        % Plot the solution
        figure(1); clf
        plot(t,y,'-*','LineWidth',2);
        legend({'y_{h1}', 'y_{h2}', 'y_{h3}', 'y_{h4}'});
        hold on
        plot(t,sol_t,'-','LineWidth',1);
        xlabel('t');
        ylabel('y_h(t)');
        title(sprintf('Solution %slicit Euler for h = %8.2e, %s grid',exp_imp,h,grid));
        %   figfontifier(gcf,15);
        count=count+1;
        if (savefigs==1)&&(count==2)
            saveas(gcf,sprintf('demo_stiff_sol_%s_%s',exp_imp,grid),'epsc');
            axis([0 T 1.5*min(min(sol_t)) 1.5*max(max(sol_t))]);
            saveas(gcf,sprintf('demo_stiff_sol_%s_%s_clipp',exp_imp,grid),'epsc');
            if strcmp(grid,'switched')
                axis([0 5*ts 1.5*min(min(sol_t)) 1.5*max(max(sol_t))]);
                saveas(gcf,sprintf('demo_stiff_sol_%s_%s_zoom',exp_imp,grid),'epsc');
            end
        end
        
        % Plot the error
        figure(2);
        maxE=1e5;
        if E<maxE
            loglog(h,E,'o','MarkerFaceColor','blue');
        else
            
            loglog(h,maxE,'^','MarkerFaceColor','red');
        end
        hold on
        axis([1/(2*max(nlist)) 1/(2*min(nlist)) 1e-5 maxE*1.5]);
        xlabel('h');
        ylabel('maximum error');
        title(sprintf('Maximum error by %slicit Euler, %s grid',exp_imp, grid));
        slope_triangle(1e-3,5e1,1e-4,1,'blue');
        %   figfontifier(gcf,15);
        
        % Wait a little
        %   pause(0.5);
    end
    
    
    % % Finally, store the error plot
    % if savefigs==1
    %   figure(2);
    % %   figfontifier(gcf,15);
    %   saveas(gcf,sprintf('demo_stiff_%s_%s',exp_imp,grid),'epsc');
    %   saveas(gcf,sprintf('demo_stiff_%s_%s',exp_imp,grid),'png');
    % end
    
end % for themethod