solve_pde


function solve_pde
    % Define the number of elements and grid points
    N = 20;  % Number of intervals
    x = [0,sort(rand(1, N-1)),1];  % Uniform grid points
    
    % Define the source term f(x)
    f = @(x) 1;  % Constant source term
    
    % Compute the finite element matrix A
    A = fem_matrix_general_grid(x);
    
    % Compute the right-hand side vector b
    b = fem_rhs_general_grid(x, f);
    
    % Apply Dirichlet boundary conditions (u(0) = u(1) = 0)
    % Remove the first and last rows/columns of A and b
    A = A(2:end-1, 2:end-1);
    b = b(2:end-1);
    
    % Solve the linear system
    u_interior = A \ b;
    
    % Reconstruct the full solution including boundary values
    u = [0; u_interior; 0];

    % Analytic solution
    c = [1, 1; exp(1), exp(-1)]\[-1; -1];
    u_analytic = @(x) c(1).*exp(x) + c(2).*exp(-x) + 1;
    
    % Plot the solution
    plot(x, u, '-o');
    hold on;
    plot(linspace(0,1,100), u_analytic(linspace(0,1,100)));
    xlabel('x');
    ylabel('u(x)');
    title('Solution to -u'''' + u = 1 in H^1_0((0,1))');
    grid on;
end

function A = fem_matrix_general_grid(x)
    % fem_matrix_general_grid computes the finite element matrix for
    % p(x) = 1 and q(x) = 1 on a general grid.
    %
    % Input:
    % x - a vector of grid points [x0, x1, ..., xN]
    %
    % Output:
    % A - the finite element matrix of size (N+1) x (N+1)
    
    N = length(x) - 1;  % Number of elements (intervals)
    A = zeros(N+1, N+1);  % Initialize the matrix
    
    % Loop through each element (each interval between grid points)
    for i = 1:N
        h = x(i+1) - x(i);  % Element size
        
        % Local stiffness matrix (for p(x) = 1)
        K_local = [1, -1; -1, 1] * (1/h);
        
        % Local mass matrix (for q(x) = 1)
        M_local = [2, 1; 1, 2] * (h/6);
        
        % Local matrix is the sum of stiffness and mass contributions
        A_local = K_local + M_local;
        
        % Add local matrix to the global matrix
        A(i:i+1, i:i+1) = A(i:i+1, i:i+1) + A_local;
    end
end

function b = fem_rhs_general_grid(x, f)
    % fem_rhs_general_grid computes the right-hand side vector for
    % a variational problem with the source term f(x) on a general grid,
    % using MATLAB's built-in integration function.
    %
    % Input:
    % x - a vector of grid points [x0, x1, ..., xN]
    % f - a function handle for the source term f(x)
    %
    % Output:
    % b - the right-hand side vector of size (N+1)

    N = length(x) - 1;  % Number of elements (intervals)
    b = zeros(N+1, 1);  % Initialize the RHS vector

    % Loop through each element
    for i = 1:N
        h = x(i+1) - x(i);  % Element size

        % Define the local basis functions on the current element [x(i), x(i+1)]
        phi1 = @(xi) (x(i+1) - xi) / h;  % Hat function phi_i (left node)
        phi2 = @(xi) (xi - x(i)) / h;    % Hat function phi_{i+1} (right node)

        % Compute the local right-hand side contributions using MATLAB's integral function
        % Integral for phi1 (basis function corresponding to node i)
        b(i) = b(i) + integral(@(xi) f(xi) .* phi1(xi), x(i), x(i+1));

        % Integral for phi2 (basis function corresponding to node i+1)
        b(i+1) = b(i+1) + integral(@(xi) f(xi) .* phi2(xi), x(i), x(i+1));
    end
end