diff --git a/CDM.m b/CDM.m
new file mode 100644
index 0000000..58a1c6e
--- /dev/null
+++ b/CDM.m
@@ -0,0 +1,59 @@
+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
diff --git a/Plot1.PNG b/Plot1.PNG
new file mode 100644
index 0000000..8881490
Binary files /dev/null and b/Plot1.PNG differ
diff --git a/Plot2.PNG b/Plot2.PNG
new file mode 100644
index 0000000..2f21c46
Binary files /dev/null and b/Plot2.PNG differ
diff --git a/defl.asv b/defl.asv
new file mode 100644
index 0000000..365d0c6
--- /dev/null
+++ b/defl.asv
@@ -0,0 +1,23 @@
+function dwdt = defl(x,w)
+%For ODE45 script - recursively defines fourth order ODE for deflection
+%given q=0 and P=0
+
+%Initialize Constants
+E = 70*10^9; %Pa
+dens = 2700; %kg/m^3
+b = 0.1; %m
+h = 0.01; %m
+l = 1; %m
+
+%Placeholder frequenc
+%freq = 150;
+
+%Derived variables
+I = (b*(h^3))/12; % m^4
+area = b*h; %m^2
+
+%Define constant K in place
+k = (dens*area*(freq^2))/(E*I);
+dwdt = [w(2);w(3);w(4);k*w(1)];
+end
+
diff --git a/defl.m b/defl.m
new file mode 100644
index 0000000..a73df33
--- /dev/null
+++ b/defl.m
@@ -0,0 +1,24 @@
+function dwdt = defl(x,w)
+%For ODE45 script - recursively defines fourth order ODE for deflection
+%given q=0 and P=0, x and w used by ODE45 to solve for deflection over
+%position
+
+%Initialize Constants
+E = 70*10^9; %Pa
+dens = 2700; %kg/m^3
+b = 0.1; %m
+h = 0.01; %m
+l = 1; %m
+
+%Placeholder frequency
+%freq = 150;
+
+%Derived variables
+I = (b*(h^3))/12; % m^4
+area = b*h; %m^2
+
+%Define constant K to simplify recursive functions for ODE
+k = (dens*area*(freq^2))/(E*I);
+dwdt = [w(2);w(3);w(4);k*w(1)];
+end
+
diff --git a/egv.m b/egv.m
new file mode 100644
index 0000000..c432cbc
--- /dev/null
+++ b/egv.m
@@ -0,0 +1,61 @@
+function [f1,f2,f3] = egv(N,P)
+%Function solves for the natural frequency of a the beam in problem 2 and 3
+%given N segments and transverse load P
+
+%Initialize Constants
+E = 70*10^9; %Pa
+dens = 2700; %kg/m^3
+b = 0.1; %m
+h = 0.01; %m
+l = 1; %m
+
+%Derived variables
+I = (b*(h^3))/12; % m^4
+area = b*h; %m^2
+
+%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
+
+%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;
+
+%Built in eigenvalue MATLAB function finds eigenvalues of the matrix
+%generated for the beam
+egv = eig(A);
+
+%Plug into function for natural frequencies to solve (for first three -
+%more can be generated with following values of egv vector)
+f1 = sqrt((egv(1)*E*I)/dens/area/(segl^4));
+f2 = sqrt((egv(2)*E*I)/dens/area/(segl^4));
+f3 = sqrt((egv(3)*E*I)/dens/area/(segl^4));
+end
diff --git a/montecarlo.m b/montecarlo.m
new file mode 100644
index 0000000..2e1bfe6
--- /dev/null
+++ b/montecarlo.m
@@ -0,0 +1,30 @@
+function [avg,stdev]=montecarlo(num)
+%Generates random values for base and height values
+%of the beam, finds according mean and standard deviation
+%of the outputted deflections.
+%Applied distributed load given as 50 N/m.
+%Normal random number generation using 0.1% standard deviation.
+%Input 'num' designates number of random samples.
+
+%Given physical parameters
+q = 50; %Input loading 50 N/m
+E = 70*10^9; %Pa 
+l = 1; %m
+x = 0.5; %m
+
+%Loop generates random numbers and places calculated deflections in matrix
+for i=1:num
+    %Generates random values for I
+    b = normrnd(0.1,0.01); %random base length
+    h = normrnd(0.01,0.001); %random height value
+    I = (b*h^3)/12; %calculates random Moment of Inertia (m^4)
+
+    %Calculated corresponding deflections and place in matrix w
+    w(i) = -(((q*l*(x^3))/(12*E*I)))-(((q*(x^4))/(24*E*I)))-((q*(l^3)*x)/(24*E*I));
+end
+
+%Find descriptive statistics for set of delfections
+avg = mean(w);
+stdev = std(w);
+end
+