%% 3 Finite element method with Gauss quadrature

figure();
hold('on');

for n=[10,20,40,80,160]
    L=10;
    h=L/n;
%why is this this way? (the two last lines are extra)
    A=2/h*diag(ones(n-1,1),0)...
        -1/h*diag(ones(n-2,1),-1)...
        -1/h*diag(ones(n-2,1),1)...
        -1/2*diag(ones(n-2,1),-1)...
        +1/2*diag(ones(n-2,1),1);

    b=zeros(n-1,1);
    f=@(x) 5*exp(x);
    phi_l=@(z) (1+z)/2;
    phi_r=@(z) (1-z)/2;
    for i=1:n-1
        z2x_l=@(z) (i-1/2)*h+z*h/2;
        g_l=@(z) phi_l(z)*f(z2x_l(z));
        z2x_r=@(z) (i+1/2)*h+z*h/2;
        g_r=@(z) phi_r(z)*f(z2x_r(z));
        b(i)=h/2*(1*g_l(1/sqrt(3))+1*g_l(-1/sqrt(3)))+h/2*(1*g_r(1/sqrt(3))+1*g_r(-1/sqrt(3)));
    end

    u=[0;A\b;0];
    x=0:h:L;
    plot(x,u,'k-');
end

ua=-exp(x).*(5*x-(5*(1+9*exp(10)))/(-1+exp(10))-5)-50*exp(10)/(-1+exp(10));
plot(x,ua,'r-');

hold('off');
xlabel('x');
ylabel('u');