CDM.m

function [w] = CDM(q,N,P)
%Given an applied distributed load (q N/m), N segments, and P newtons
%transverse load, uses Central Difference Method to solve for deflection at
%segments

%Initialize Constants
E = 70*10^9; %Pa
b = 0.1; %m
h = 0.01; %m
l = 1; %m

%Derived variables
I = (b*(h^3))/12; % m^4

%Segments
segl = l/N; %m - length of segment

%Diagonal values
diagonal = ((P*2)/((N^2)*E*I)) + 6; %main diagonal
offDiag = (-P/((N^2)*E*I)) - 4; %off-diagonal

%Deflection differential 
dw = ((segl^4)*q)/(E*I);

%Loop through all values of A to create values
A = zeros(N-1,N-1);
for row = 1:(N-1)
    for column = 1:(N-1)
        %diagonal values
        if row == column
            A(row,column) = diagonal;
        end
        %off diagonal values
        if column == row - 1
            A(row,column) = offDiag;
        end
        if column == row + 1
            A(row,column) = offDiag;
        end
        %"off-off" diagonal values
        if column == row - 2
            A(row,column) = 1;
        end
        if column == row + 2
            A(row,column) = 1;
        end
    end
end

%Sets corner values to one less than main diagonal values
A(1,1) = ( (2*P) / ((N^2)*(E*I)) ) + 5;
A(N-1,N-1) = ( (2*P) / ((N^2)*(E*I)) ) + 5;

%Generate B Matrix
B = dw*ones((N-1),1);

% Calculate deflection
w = A\B;
end