%Question1
clc;clear all; close all;
syms x u1 u2 u3 u4 u5 Phi I;
Phi(1,1) = 1-4*x;
Phi(2,1) = 4*x; Phi(2,2) = 2-4*x;
Phi(3,2) = 4*x-1; Phi(3,3) = 3-4*x;
Phi(4,3) = 4*x-2; Phi(4,4) = 4-4*x;
Phi(5,4) = 4*x-3;

coord = [0, 0.25, 0.5, 0.75, 1];
for i=1:5 % Integral loop
    I(i) = 0; % Initialize the i-th integral to zero.
    for j=1:4 % Element loop
        w = Phi(i,j);
        dwdx = diff(w, x);
        u = u1*Phi(1,j) + u2*Phi(2,j) + u3*Phi(3,j) + u4*Phi(4,j) + u5*Phi(5,j);
        dudx = diff(u, x);
        I(i) = I(i) + int(dwdx*dudx - w*u + w*x^2, x, coord(j), coord(j+1));
    end
end

% Linear equation system matrix
I = [
    (47/6),    (-97/24),    0;
    (-97/24),  (47/6),    (-97/24);
    0,         (-97/24),  (47/6)
    ];

% Define the right-hand side vector
b = [-7/384 ; -25/384 ;-55/384];

% Solve the linear equation system
u = I\b

% Plotting
coord = [0, 0.25, 0.5, 0.75, 1];
u = [0 -0.0232 -0.0405 -0.0392 0]; % Assuming these are the correct values of u
plot(coord, u, 'o-b', 'LineWidth', 2)
hold on

x = 0:0.01:1;
uExact = (sin(x) + 2*sin(1-x)) / sin(1) + x.*x - 2;
plot(x, uExact, 'r', 'LineWidth', 2)

%-------------------------------------------------------------