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);
z = 1:10;

for k = 1:10
    n = (10 * 2^k)-1; %+1 so that x=L/2 is included in our domain 
    g = 9.81; %downward force due to gravity
    h = L/(n+1); %define variables
    midpoint = (n/2)+.5; %get midpoint of domain
    x = h:h:L-h;
    I = w*d^3/12;
    E = 1.3e10;
    p = 100;

A = zeros(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) = [1 -4 6 -4];
A(n,n-3:n) = [-1/4 8/3 -9 16];

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-2*L.*x.^3+L^2.*x.^2)-((p*g*L^2)/(pi^4*E*I).*(L^2*sin(pi/L.*x)+pi.*x.*(x-L))))';
    error = abs(ysolution - yactual);
    errors(k) = error(midpoint);
end


errors
figure();
semilogy (nvalue, errors, 'ro', nvalue, errors);
title('Error')
xlabel('n');
ylabel('error(n)');

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

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

table(nvalue,errors)

table(nvalue,condition)