function MainFile

clc

N = 10; alpha = 0.1;
x = zeros(N,1);
A1 = zeros(N,1); B1 = zeros(N,1); R1 = zeros(N,1);
A2 = zeros(N,1); B2 = zeros(N,1); R2 = zeros(N,1);
R = zeros(N,1); c = ones(N,1); J = zeros(N,N);

for k = 1:N;
    x(k+1) = x(k)+(1/N);
end

% disp(x)

for i = 2 %looping over elements
    A1(i) = c(i)*2*(0.25*(0.5*x(i+1)-0.5*x(i))^-1)+c(i+1)*2*(-0.25*(0.5*x(i+1)-0.5*x(i))^-1);
    B1(i) = (alpha*0.5*(x(i+1)-x(i))*0.125*(c(i))^2)*((5/9)*(1+sqrt(3/5))^3+(8/9)+(5/9)*(1-sqrt(3/5)))+...
        (alpha*0.5*(x(i+1)-x(i))*0.25*c(i)*c(i+1))*((5/9)*(1-sqrt(3/5)*(1+sqrt(3/5))^2)+(8/9)+(5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2)+...
        (alpha*0.125*c(i+1)^2*0.5*(x(i+1)-x(i)))*((5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2+(8/9)+(5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2);
    R1(i) = A1(i)+B1(i);
    
    A2(i) = c(i)*2*(-0.25*(0.5*x(i+1)-0.5*x(i))^-1)+c(i+1)*2*(0.25*(0.5*x(i+1)-0.5*x(i))^-1);    
    B2(i) = (alpha*0.5*(x(i+1)-x(i))*0.125*(c(i))^2)*((5/9)*(1-sqrt(3/5))^3+(8/9)+(5/9)*(1+sqrt(3/5)))+...
        (alpha*0.5*(x(i+1)-x(i))*0.25*c(i)*c(i+1))*((5/9)*(1+sqrt(3/5)*(1-sqrt(3/5))^2)+(8/9)+(5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2)+...
        (alpha*0.125*c(i+1)^2*0.5*(x(i+1)-x(i)))*((5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2+(8/9)+(5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2);
    R2(i) = A2(i)+B2(i);
    
    R(i) = R1(i) + R2(i);
    
    J(1,1) = 1;
   
    J(i,i) = 2*(0.25*(0.5*x(i+1)-0.5*x(i))^-1)+c(i+1)*2*(-0.25*(0.5*x(i+1)-0.5*x(i))^-1)+...
        (2*alpha*0.5*(x(i+1)-x(i))*0.125*(c(i)))*((5/9)*(1+sqrt(3/5))^3+(8/9)+(5/9)*(1-sqrt(3/5)))+...
        (alpha*0.5*(x(i+1)-x(i))*0.25*c(i+1))*((5/9)*(1-sqrt(3/5)*(1+sqrt(3/5))^2)+(8/9)+(5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2)+...
        (alpha*0.125*c(i+1)^2*0.5*(x(i+1)-x(i)))*((5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2+(8/9)+(5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2);
    
    J(i,i+1) = c(i)*2*(0.25*(0.5*x(i+1)-0.5*x(i))^-1)+2*(-0.25*(0.5*x(i+1)-0.5*x(i))^-1)+...
        (alpha*0.5*(x(i+1)-x(i))*0.125*(c(i))^2)*((5/9)*(1+sqrt(3/5))^3+(8/9)+(5/9)*(1-sqrt(3/5)))+...
        (alpha*0.5*(x(i+1)-x(i))*0.25*c(i))*((5/9)*(1-sqrt(3/5)*(1+sqrt(3/5))^2)+(8/9)+(5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2)+...
        (2*alpha*0.125*c(i+1)*0.5*(x(i+1)-x(i)))*((5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2+(8/9)+(5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2);

    J(i+1,i) = 2*(-0.25*(0.5*x(i+1)-0.5*x(i))^-1)+c(i+1)*2*(0.25*(0.5*x(i+1)-0.5*x(i))^-1)+...
        (2*alpha*0.5*(x(i+1)-x(i))*0.125*(c(i)))*((5/9)*(1-sqrt(3/5))^3+(8/9)+(5/9)*(1+sqrt(3/5)))+...
        (alpha*0.5*(x(i+1)-x(i))*0.25*c(i+1))*((5/9)*(1+sqrt(3/5)*(1-sqrt(3/5))^2)+(8/9)+(5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2)+...
        (alpha*0.125*c(i+1)^2*0.5*(x(i+1)-x(i)))*((5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2+(8/9)+(5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2);

    J(i+1,i+1) = c(i)*2*(-0.25*(0.5*x(i+1)-0.5*x(i))^-1)+2*(0.25*(0.5*x(i+1)-0.5*x(i))^-1)+...
        (alpha*0.5*(x(i+1)-x(i))*0.125*(c(i))^2)*((5/9)*(1-sqrt(3/5))^3+(8/9)+(5/9)*(1+sqrt(3/5)))+...
        (alpha*0.5*(x(i+1)-x(i))*0.25*c(i))*((5/9)*(1+sqrt(3/5)*(1-sqrt(3/5))^2)+(8/9)+(5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2)+...
        (2*alpha*0.125*c(i+1)*0.5*(x(i+1)-x(i)))*((5/9)*(1-sqrt(3/5))*(1+sqrt(3/5))^2+(8/9)+(5/9)*(1+sqrt(3/5))*(1-sqrt(3/5))^2);
        
end

% disp(R1)
% disp(R2)
% disp(J)
disp(R)

% 
% for COUNTER = 1:100
%     delta_c = J\(-R);
%     c = c+delta_c;
%     
%     normRatio = norm(delta_c)/norm(c);
%     
%     if (normRatio)< 1E-5
%         disp ('Newtons works');
%         break
%     end
%     
%     if COUNTER ==100
%         disp('Newton FAIL');
%         break
%     end
% end
% 
% end