% Code for implementing Galerkin Method to the general 1-D 2nd order ODE;
% the problem discussed in chapter 6

clc; clear all;

%-------Input-----------------%
n = input('Number of elements: ');

%-------Calculation---------%
L = 1;
Le = L/n;      %length of one element

x = ones(n+1, 1);            %to store x-coordinates


k_e = zeros(2,2);             %Local 2 x 2 stiffness matrix
f_e = zeros(2,1);             %Local 2 x 1 force vector
K_G = zeros(n+1,n+1);         %Global stiffness matrix of size DOF x DOF
F_G = zeros(n+1, 1);          %Global force matrix

for i = 1:n
    x_1 = 1+(i-1)*Le;                 % cordinate of left of the element
    x_2 = x_1+Le;                     % cordinate of right of the element
    x(i+1) = x_2;                     %x for node 2
        
    k(i) = 0.5*(x_2^2-x_1^2)/(x_2-x_1)^2;           %spring constant of the bar element
    k_e = k(i)*[1 -1;-1 1];                         %local stiffness matrix
    
    f_e(1) = -4/(x_2-x_1)*(x_2^3/6-x_2*x_1^2/2+x_1^3/3);  %constant part of f_e matrix, the derivative is left for calculation
    f_e(2) = -4/(x_2-x_1)*(x_1^3/6-x_1*x_2^2/2+x_2^3/3);
    
    %------Assembly to global matrix----------%
   K_e = zeros(n+1, n+1);        %Local stiffness matrix with additional rows and columns of zeros
   K_e([i:i+1],[i:i+1])=k_e;
   K_G = K_G+K_e; 
   F_G(i)=F_G(i)+f_e(1);
   F_G(i+1)=F_G(i+1)+f_e(2);
end

%------------------Solution-----------------------%


K_G_act = K_G([2:n],[2:n]);    %Since node 1 and the last node's displacements are zero
F_G_act = F_G([2:n]);          %First and the last equations will be the part of the constraint equation
Y_act = inv(K_G_act)*F_G_act;      %Calculation of unknown displacements
Y_tot = [0; Y_act;0];                %all displacements

plot(x,Y_tot);                   
xlabel('Y');
ylabel('X');