clear; clc;
lambda = @(x) (x < 1/sqrt(2)) + 10*(x >= 1/sqrt(2)); %diffusion coefficient 
f = @(x) 0; %homogeneous RHS
x = linspace(0,1,6);   %coarse initial mesh, doesn't include x_m
nIter = 6;

for it = 1:nIter
    [K,F] = assemble_1D(x,lambda,f); %assemble
    K(1,:) = 0;  K(1,1) = 1;  F(1) = 0; % Dirichlet BC at x=0
    K(end,:) = 0; K(end,end) = 1; F(end) = 1; % Dirichlet BC at x=1
    u = K\F;

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

figure
plot(x,u,'-o','LineWidth',1.5)
xlabel('x'), ylabel('u')
title('Ex6B, h-adaptivity with coefficient jump')
grid on