centraldifference2.m

function centraldifference2(n)
%Takes the number of segments as the input argument and then runs a for 
%loop of the P value and a for loop of the q value to find the maximum deflection
%with those given parameters, the calculated max deflcetion valus are then
%plotted with respect to the given distributed loads.
%Within the for loop the matrix is solved for using the central difference
%formulas and patterns that it follows and then the backslash operator is
%used
b=0.1;
h=0.01;
I=b*(h^3)/12;
E=70e9; 

P=[0,100,200,300];
q=[1,10,20,30,50];
for p=1:4 
    for i=1:5
        D=-4-((P(p)*((1/n)^2))/(E*I));
        e=ones((n-2),1)*D;
        F= 6+(P(p)*(2*(1/n)^2)/(E*I));
        f=ones(n-1,1)*F;
        g=ones((n-3),1);
        A=diag(e,-1)+ diag(f)+diag(e,1)+diag(g,2)+diag(g,-2);
        A(1,1)=5+(P(p)*(2*(1/n)^2)/(E*I));
        A((n-1),(n-1))=5+(P(p)*(2*(1/n)^2)/(E*I));   
        b=(((1/n)^4)*q(i))/(E*I);
    if n==6
    B=b*ones((n-1),1); 
    W=A\B;
    dx(i)=max(W);
elseif n==10
    B=b*ones((n-1),1);
    W=A\B;
    dx(i)=max(W);
elseif n==20
    B=b*ones((n-1),1);
    W=A\B;
    dx(i)=max(W);
    end
   
    end
    DX(p,:)=dx;
 plot(q,DX)
xlabel('q (force/length)')
ylabel('dx (max deflection)')
title('Load vs Max Deflection')
legend('P=0','P=100','P=200','P=300')
end
end