clc
clear all
close all


% coordinates for visualization
xx = linspace(0,1,100)';

% analytical solution at points xx:
u =  -2*xx.^2 + 4*xx - 4;

% number of nodes in different meshes
Nodes = [4 6 10 20 40 50];

% solve the problem with different number of nodes:
for ii = 1:length(Nodes)

    % collect number of nodes
    N = Nodes(ii);

    % mesh size
    h = 1/(N-1);

    % this is the FE grid

    % x-positions of nodes
    g  = (0:h:1)';
    h1 = (1:N-1)';  h2 = (2:N)';

    % element topology
    H  = [h1 h2];
    clear h1 h2

    % solve the FE-equation

    % matrices
    S = zeros(N,N);
    M = zeros(N,N);

    % loop over elements
    for kk = 1:size(H,1)

        % collect node indices for element kk
        nodes = H(kk,:);

        % vertices for element kk
        xkk = g(nodes,:);

        % h for element kk (constant for each element but does not need to be)
        hkk = xkk(2) - xkk(1);

        % S matrix for linear element
        S_kk = [1 -1;-1 1] / hkk;

        % M matrix for a linear element
        M_kk = [1/3, 1/6; 1/6, 1/3] * hkk;

        % collect
        S(nodes,nodes) = S(nodes,nodes) + S_kk;
        M(nodes,nodes) = M(nodes,nodes) + M_kk;
    end
    clear M_kk S_kk

    % collect meshsize (constant)
    h_ell(ii) = hkk;

    % visualize matrices
    figure(1)
    subplot(1,2,1)
    spy(S)
    title('S matrix')
    subplot(1,2,2)
    spy(M)
    title('M matrix')


    % form system matrix
    K = S + M;
    clear S

    % bc term (phi_j(0) = 1, otherwise this is zero)
    bc    = zeros(N,1);
    bc(1) = 1;

    % original right hand side time (basis functions are ones on those node
    % points)
    f = -2*g.^2 + 4*g;

    % right hand side
    b = M*f - 4*bc;
    clear M

    % estimate using backslash
    alphavec = K\b;

    % calculate the solution at points xx:
    uh = zeros(length(xx),1);
    for jj = 1:length(xx)
        xxjj = xx(jj);
        Eljj = [];
        nodejj = find(g==xxjj);
        if ~isempty(nodejj) % point tjj is a node
            uh_jj = alphavec(nodejj);
        else                 % point tjj is not a nodal point
            for kk = 1:size(H,1)  % find the element Eljj which includes point xxjj
                nodes = H(kk,:);
                xkk = g(nodes,:);
                if xkk(1) < xxjj & xxjj < xkk(2)
                    Eljj = kk;
                    break
                end
            end

            % let's find coefficients those satisfy:
            % y(x = x_0) = a * x + b = alpha_0
            % y(x = x_1) = a * x + b = alpha_1
            locMat = [xkk ones(2,1)];
            loc_coeffs = locMat \ alphavec(nodes);
            % now we can evaluate the basis inside the element
            uh_jj = loc_coeffs(1)*xxjj + loc_coeffs(2);

        end

        % collect value
        uh(jj) = uh_jj;
    end

    % calculate the l2 error:
    deltauh_ii = sqrt(sum((uh-u).^2));
    deltauh(ii) = deltauh_ii;

    % relative l2-error
    u_ii = sqrt(sum((u).^2));
    relerr(ii) = deltauh_ii/u_ii;

    % analytical solution and FE-solution
    figure(2)
    subplot(2,3,ii)
    hold on
    p1 = plot(xx,u,'k');
    p3 = plot(xx,uh,'r');
    set(p1,'LineWidth',2);

    xlabel('$x$', 'Interpreter','latex','FontSize',12);
    ylabel('solution amplitude', 'Interpreter','latex','FontSize',12);
    grid on
    xlim([0 1])
    ylim([-4.5 -2])

    hold off
    title(['N = ',num2str(N)])
    pause
end

% L2-error against meshsize
figure(3), clf
subplot(1,2,1)
semilogy(1./h_ell,relerr,'*-r')
axis([0 51 -0.001 0.18])
xl = xlabel('$1/h$', 'Interpreter','latex','FontSize',12);
yl = ylabel('$\|u-u_{exact}\|_2/\|u_{exact}\|_2$', 'Interpreter','latex','FontSize',12);
grid on
axis tight

subplot(1,2,2)
for ii = 1:length(Nodes)-1
    O(ii) = log((deltauh(ii+1))/(deltauh(ii)))/log(h_ell(ii+1)/h_ell(ii));
end
loglog(1./h_ell, sqrt(deltauh))
xl = xlabel('$1/h$', 'Interpreter','latex','FontSize',12);
yl = ylabel('$\|u-u_{exact}\|_2$', 'Interpreter','latex','FontSize',12);
grid on

disp('numerical convergence:')
O