% Sheet 4 / Task C
clf

nvec = [10, 20, 30, 40, 70];    % vector with number of elements
p = 70; % parameter of ode

for k=1:length(nvec)
    ne = nvec(k);
    nodes = linspace(0,1,ne+1);

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

    for e = 1:ne
        xleft = nodes(e);
        xright = nodes(e+1);

        % first part of the bilinearform phi'*phi'
        funcone = @(x) 1 +0*x;
        Ke = ne^2*[integral(funcone,xleft,xright), -integral(funcone,xleft,xright);
            -integral(funcone,xleft,xright), integral(funcone,xleft,xright)];

        % second part: phi'*phi
        func1 = @(x) xright - x;
        func2 = @(x) x - xleft;
        Ke = Ke + p*ne^2*[-integral(func1,xleft,xright), integral(func1,xleft,xright);
                        -integral(func2,xleft,xright), integral(func2,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;

    hold on
    plot(nodes,u)
    title('approximation of solution u for p=70');
    xlabel('x values');
    ylabel('u_h');
    grid on;
    leg = legend(num2str(nvec'),'Location', 'Northwest');
    title(leg, 'n =');
     
end
saveas(gcf, 'ex_4_C_plot.jpg');