% Sheet 4 / Task B

ne = 50;    % number of elements
nodes = linspace(0,1,ne+1);

K = zeros(ne+1,ne+1);   % global stiffness matrix
f = zeros(ne+1,1);    % load vector

lambda = @(x) (x<1/sqrt(2))*1 + (x>=1/sqrt(2))*10;

for e = 1:ne
    xleft = nodes(e);
    xright = nodes(e+1);
    Ke = ne^2*[integral(lambda,xleft,xright), -integral(lambda,xleft,xright);
        -integral(lambda,xleft,xright), integral(lambda,xleft,xright)];
    
    K(e:e+1, e:e+1) = K(e:e+1, e:e+1) + Ke;
end

% disp(K)

% adaption for dirichlet boundary
K(1,1) = K(1,1)*(1 + 1e6);
f(end) = K(end,end)*1e6;
K(end,end) = K(end,end)*(1 + 1e6);

u = K\f;

plot(nodes,u)
title('approximation of solution u');
xlabel('x values');
ylabel('u_h');