clear; clc; close all;
pList = [5, 10, 20, 100];
nIter = 5;

figure;
tiledlayout(2,2)
for ip = 1:length(pList)
    p = pList(ip);

    f = @(x) 2*p^3*x./(p^2*x.^2+1).^2; %RHS
    lambda = @(x) 1; %DIffusion coefficient 

    u_exact = @(x) atan(p*x);

    x = linspace(-1,1,10); %initial uniform mesh (includes x=0)

    for it = 1:nIter %adaptive FEM loop
        [K,F] = assemble_1D(x,lambda,f); %assemble global stiffness matrix K and load vector F
        %Neumann boundary conditions at x=1
        F(end) = F(end) + p/(p^2+1);
        %Dirichlet boundary conditions at x=-1
        K(1,:) = 0;
        K(1,1) = 1;
        F(1) = -atan(p);
        % GH: Dirichlet boundary conditions at x=+1
        K(end,:) = 0;
        K(end,end) = 1;
        F(end) = atan(p);

        u = K \ F;

        if it < nIter % h-adaptivity loop
            eta = flux_jump(x,u,lambda); 
            x = h_adapt(x,eta,0.3);
        end
    end

    u_ex = u_exact(x(:));        % force column vector
    err_inf = max(abs(u - u_ex));
    nexttile
    plot(x,u,'-o','LineWidth',1.5); hold on
    plot(x,u_ex,'--','LineWidth',1.5)
    title(['p = ',num2str(p)])
    xlabel('x'), ylabel('u')
    legend('FEM','exact','Location','best')
    grid on

    fprintf('p = %d\n', p);
    fprintf('error = %.3e\n\n', err_inf);
end

sgtitle('Exercise 6A: h-adaptivity and exact solution comparison')