%% --- Complete Analysis of Implicit Euler for y' + Ay = 0 ---
clear; close all; clc;

%% 1. Problem Definition
A = [2 1; 1 3];      % 2x2 SPD matrix for phase portraits
y0 = [2; 3];         % Initial condition
n = size(A, 1);      % Dimension of the system

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2. Convergence Test (Your original code)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('--- Running Convergence Test ---\n');

% --- Parameters for Convergence Test ---
T_conv = 1.0;                     % Final time for convergence test
y_exact_conv = expm(-A*T_conv) * y0; % Exact solution at T_conv

% --- Run Test for Different Step Sizes ---
Nlist = [10 20 40 80 160 320];
hlist = T_conv ./ Nlist;
errors = zeros(size(Nlist));

for k = 1:length(Nlist)
    N = Nlist(k);
    h = hlist(k);
    
    % Initialize solution at y0
    y = y0;
    
    % Define the implicit update matrix
    M = (eye(n) + h*A);
    
    % Time-stepping loop
    for j = 1:N
        y = M \ y; % Implicit Euler update
    end
    
    % Calculate the L2 norm of the error at the final time
    errors(k) = norm(y - y_exact_conv);
end

% --- Plot Convergence Results ---
figure;
loglog(hlist, errors, '-o', 'LineWidth', 1.5, 'MarkerFaceColor', 'b');
hold on;
loglog(hlist, hlist, 'k--', 'LineWidth', 1); % Plot a reference line with slope 1
xlabel('Step size (h)');
ylabel('Error at T=1');
title('Convergence of Implicit Euler');
legend('Error', 'Reference (Order 1)', 'Location', 'northwest');
grid on;
axis tight;

% --- Display Convergence Table ---
fprintf('Convergence table:\n');
conv_table = table(Nlist', hlist', errors', 'VariableNames', {'N','h','Error'});
disp(conv_table);
fprintf('Note: The slope on the log-log plot is ~1, confirming first-order accuracy.\n\n');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3. Stability Demonstration (Using Phase Portraits)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('--- Running Stability Demonstration ---\n');

% --- Parameters for Stability Test ---
T_stab = 20;         % Long final time to observe stability

% --- Determine Stability Limit for Explicit Euler ---
lambda = eig(A);
lambda_max = max(lambda);
h_limit_explicit = 2 / lambda_max;
fprintf('Explicit Euler is stable only if h < %.3f\n', h_limit_explicit);

% --- Choose a step size LARGER than this limit ---
h_stab = h_limit_explicit * 1.1; 
fprintf('Using h = %.3f for the demonstration.\n', h_stab);

% --- Simulate with Both Methods ---
N_stab = round(T_stab / h_stab);
t_stab = linspace(0, T_stab, N_stab+1);
y_hist_implicit = zeros(n, N_stab+1);
y_hist_explicit = zeros(n, N_stab+1);
y_hist_exact    = zeros(n, N_stab+1);

y_imp = y0; y_exp = y0;
y_hist_implicit(:, 1) = y_imp;
y_hist_explicit(:, 1) = y_exp;

M_implicit = (eye(n) + h_stab*A);
M_explicit = (eye(n) - h_stab*A);

for j = 1:N_stab
    y_imp = M_implicit \ y_imp;
    y_exp = M_explicit * y_exp;
    y_hist_implicit(:, j+1) = y_imp;
    y_hist_explicit(:, j+1) = y_exp;
end

for j = 1:N_stab+1
    y_hist_exact(:, j) = expm(-A*t_stab(j)) * y0;
end

% --- Plot Global View (Shows divergence) ---
figure;
hold on;
plot(y_hist_exact(1,:), y_hist_exact(2,:), 'k--', 'LineWidth', 2);
plot(y_hist_implicit(1,:), y_hist_implicit(2,:), 'b-o', 'LineWidth', 1.5);
plot(y_hist_explicit(1,:), y_hist_explicit(2,:), 'r-x', 'LineWidth', 1.5);
plot(y0(1), y0(2), 'go', 'MarkerSize', 10, 'MarkerFaceColor', 'g');
plot(0, 0, 'kp', 'MarkerSize', 16, 'MarkerFaceColor', 'k');
title(['Global View (h = ', num2str(h_stab, 2), ')']);
xlabel('y_1'); ylabel('y_2');
legend('Exact', 'Implicit Euler (Stable)', 'Explicit Euler (Unstable)', 'Start', 'Equilibrium');
grid on; axis equal;
hold off;

% --- Plot Zoomed-in View (Shows correct stable behavior) ---
figure;
hold on;
plot(y_hist_exact(1,:), y_hist_exact(2,:), 'k--', 'LineWidth', 2);
plot(y_hist_implicit(1,:), y_hist_implicit(2,:), 'b-o', 'LineWidth', 1.5);
plot(y0(1), y0(2), 'go', 'MarkerSize', 10, 'MarkerFaceColor', 'g');
plot(0, 0, 'kp', 'MarkerSize', 16, 'MarkerFaceColor', 'k');
title(['Zoomed-in View (h = ', num2str(h_stab, 2), ')']);
xlabel('y_1'); ylabel('y_2');
legend('Exact', 'Implicit Euler (Stable)', 'Start', 'Equilibrium');
grid on; axis equal;
axis_limit = max(abs(y0)) * 1.5;
xlim([-axis_limit, axis_limit]);
ylim([-axis_limit, axis_limit]);
hold off;

fprintf('--- Analysis Complete ---\n');