1D finite element method in the Hermite basis (P3C1) - Problem of solution calculation

Hello to all!

I am currently working on solving the problem and Neumann conditions on the edge, with the finite element method.

I managed to set up a code that works for P1 and P2 Lagragne finite elements (i.e: linear and quadratic) and the results are good!

I am trying to implement the finite element method using the Hermite basis. This basis is defined by the following basis functions and derivatives:

syms x
phi(x) = [2*x^3-3*x^2+1,-2*x^3+3*x^2,x^3-2*x^2+x,x^3-x^2]

phi(x) = 

fplot(phi,[0,1])

% Derivative
dphi = [6*x.^2-6*x,-6*x.^2+6*x,3*x^2-4*x+1,3*x^2-2*x]

dphi = 

fplot(dphi,[0,1])

The problem with the following code is that the solution vector u is not good. I know that there must be a problem in the S and F element matrix calculation loop, but I can't see where even though I've been trying to make changes for a week.

Can you give me your opinion? Hopefully someone can see my error.

Thanks a lot,

% -alpha*u'' + beta*u = f
% u'(a) = bd1, u'(b) = bd2;

a = 0;
b = 1;
f = @(x) (1);

alpha = 1;
beta = 1;

% Neuamnn boundary conditions
bn1 = 1;
bn2 = 0;

syms ue(x)
DE = -alpha*diff(ue,x,2) + beta*ue == f;
du = diff(ue,x);

BC = [du(a)==bn1, du(b)==bn2];

ue = dsolve(DE, BC);

figure
fplot(ue,[a,b], 'r', 'LineWidth',2)

N = 2;
nnod = N*(2+2); % Number of nodes
neq = nnod*1; % Number of equations, one degree of freedom per node
xnod = linspace(a,b,nnod);

nodes = [(1:3:nnod-3)', (2:3:nnod-2)', (3:3:nnod-1)', (4:3:nnod)'];

phi = @(xi)[2*xi.^3-3*xi.^2+1,2*xi.^3+3*xi.^2,xi.^3-2*xi.^2+xi,xi.^3-xi.^2];
dphi = @(xi)[6*xi.^2-6*xi,-6*xi.^2+6*xi,3*xi^2-4*xi+1,3*xi^2-2*xi];

% Here, just calculate the integral using gauss quadrature..
order = 5;
[gp, gw] = gauss(order, 0, 1);

S = zeros(neq,neq);
M = S; 
F = zeros(neq,1);

for iel = 1:N
    %disp(iel)
    inod = nodes(iel,:);
    xc = xnod(inod);
    h = xc(end)-xc(1);
    Se = zeros(4,4);
    Me = Se;
    fe = zeros(4,1);
    for ig = 1:length(gp)
        xi = gp(ig);
        iw = gw(ig);
        
        Se = Se + dphi(xi)'*dphi(xi)*1/h*1*iw;
        Me = Me + phi(xi)'*phi(xi)*h*1*iw;

        x = phi(xi)*xc';
        fe = fe + phi(xi)' * f(x) * h * 1 * iw;
    end
    % Assembly
    S(inod,inod) = S(inod, inod) + Se;
    M(inod,inod) = M(inod, inod) + Me;
    F(inod) = F(inod) + fe;
end
S = alpha*S + beta*M;

g = zeros(neq,1);
g(1) = -alpha*bn1;
g(end) = alpha*bn2;

alldofs = 1:neq;

u = zeros(neq,1); %Pre-allocate

F = F + g;

u(alldofs) = S(alldofs,alldofs)\F(alldofs)

Warning: Matrix is singular to working precision.

u = 8×1

NaN NaN NaN NaN NaN NaN NaN NaN

figure
fplot(ue,[a,b], 'r', 'LineWidth',2)
hold on
plot(xnod, u, 'bo')

for iel = 1:N
    inod = nodes(iel,:);
    xc = xnod(inod);
    U = u(inod);
    xi = linspace(0,1,100)';
    
    Ue = phi(xi)*U;
    Xe = phi(xi)*xc';
    plot(Xe,Ue,'b -')
end

% Gauss function for calculate the integral
function [x, w, A] = gauss(n, a, b)
    n = 1:(n - 1);
    beta = 1 ./ sqrt(4 - 1 ./ (n .* n));
    
    J = diag(beta, 1) + diag(beta, -1); 
    [V, D] = eig(J);
    x = diag(D);

    A = b - a;

    w = V(1, :) .* V(1, :);
    w = w';
    x=x';
end