Bonus_Problem_Coefficients

function [ w1 ] = Coefficients( n,EIpa,L,P,E,I )
%Coefficients calculates the values of the w matrix where the values
%represent the coefficients of each node in the computational domain. The
%coefficient matrix is of the form: Ax = b
%The problem is posed as a Ax = (lambda)x = 0 so this makes sense to write
%the matrix as a coefficient matrix

% Change in x for central difference formula for 6 segments, dx6
dx = L/n;

% Precalculate denominator in central difference equation
dx2 = dx^2;
dx4 = dx^4;

% Precalculate coefficients' constant (see derivation)
K = (EIpa)/dx4;            % Central points coefficient

% Coefficient for Second Derivative term (new term in this equation)
SD = EIpa*(P/(E*I*dx2));

% Preallocate w matrix
w1 = zeros(n+1,n+1);

% Generate coefficient matrix for computational domain
% Boundary nodes w(1) and w(n+1) should equal zero
for i = 1:n+1
    if i == 1
        w1(i,i) = 6*K+2*SD;
        w1(i,i+1) = -4*K-SD;
        w1(i,i+2) = K;
    elseif i == 2
        w1(i,i-1) = -4*K-SD;
        w1(i,i) = 6*K+2*SD;
        w1(i,i+1) = -4*K-SD;
        w1(i,i+2) = K;
    elseif i == n
        w1(i,i-2) = K;
        w1(i,i-1) = -4*K-SD;
        w1(i,i) = 6*K+2*SD;
        w1(i,i+1) = -4*K-SD;
    elseif i == n+1
        w1(i,i-2) = K;
        w1(i,i-1) = -4*K-SD;
        w1(i,i) = 6*K+2*SD;
    else
        w1(i,i-2) = K;
        w1(i,i-1) = -4*K;
        w1(i,i) = 6*K;
        w1(i,i+1) = -4*K;
        w1(i,i+2) = K;
    end
end

%% Code Debug: NOT NECESSARY (CONFIRMED VIA EMAIL)
% w2 = w1;
% % If Boundary nodes are different:
% % Written using forward and backward difference equations
% % First Node
% i = 1;
%     w2(i,i) = -3*Kf;
%     w2(i,i+1) = 4*Kf;
%     w2(i,i+2) = -Kf;
% % Last Node
% i = n+1;
%     w2(i,i) = 3*Kf;
%     w2(i,i-1) = -4*Kf;
%     w2(i,i-2) = Kf;
%     
end