%% #5
clear
L = 2; 
w = .3; 
d = .03; 
errors = zeros(10,1); %initialize error matrix
condition = zeros(10,1); %initialize matrix for the condition number
hvalue = zeros(10,1);
nvalue = zeros(10,1);

for k = 1:10
    n = 10 * 2^k; 
    g = 9.81; %downward force due to gravity
    h = L/n; %define variables
    x = h:h:L;
    I = w*d^3/12;
    E = 1.3e10;
    p = 100;
    
    A = sparse(n,n); %initialize the matrix A
    A(1,1:4) = [16 -9 8/3 -1/4]; %enter the fixed values for the matrix
    A(2,1:4) = [-4 6 -4 1];
    A(n-1,n-3:n) = [16/17 -60/17 72/17 -28/17];
    A(n,n-3:n) = [-12/17 96/17 -156/17 72/17];
    
    for i = 3:n-2 %use for loops to fill in the remaining entries in matrix A
        for j = 1:n
            A(i,i-2) = 1;
            A(i,i-1) = -4;
            A(i,i) = 6;
            A(i,i+1) = -4;
            A(i,i+2) = 1;
        end
    end
    
    condition(k) = cond(A);
    hvalue(k) = h;
    nvalue(k)=n;
    forig = -480*w*d*g;
    f = (-480*w*d*g-p*g*sin(pi/L.*x))';
    fvector = f.*ones(n,1);
    RHS = (h^4/(E*I))*(fvector);
    y = zeros(n,1);
    ysolution = A\RHS;
    yactual = ((forig/(24*E*I)).*(x.^4-4*L.*x.^3+6*L^2.*x.^2)-(p*g*L)/(E*I*pi)*((L^3/pi^3)*sin((pi/L).*x)-(x.^3)/6+(L/2).*x.^2-(L^2/pi^2).*x))';
    error = abs(ysolution - yactual);
    errors(k) = error(n);
end

errors

figure();
plot (x,ysolution,x,yactual);
title('Computed Solution vs. Actual Solution for n = 10,240');
xlabel('x');
ylabel('y(x)');

figure();
semilogy (nvalue,condition,'ro',nvalue,condition);
title('Condition Number of A');
xlabel('n');
ylabel('log(condition)');

figure();
loglog(hvalue,errors, 'ro',hvalue,errors);
title('loglog Error Plot');
xlabel('log(h)');
ylabel('log(Error)');