%% FEM vs analytical solution for -u'' + a*u = f, mixed BC
clear; clc;
%% Parameters
n = 10;               % number of elements
a = 2;                % reaction term
alpha = 5;            % Robin BC
gB = 1;               % Robin BC value
f = 3;                % constant RHS
h = 1/n;
x = linspace(0,1,n+1);

%% Assemble FEM matrices
A = zeros(n+1,n+1);
F = zeros(n+1,1);

% Linear shape functions, 1D FEM
K_diff = (1/h)*[1,-1;-1,1];
K_reac = a*(h/6)*[2,1;1,2];

for e=1:n
    nodes = [e, e+1];
    A(nodes,nodes) = A(nodes,nodes) + K_diff + K_reac;
    F(nodes) = F(nodes) + f*(h/2)*[1;1]; % f*v integral approx
end

% Dirichlet BC u(0)=0
A(1,:) = 0; A(:,1) = 0; A(1,1) = 1;
F(1) = 0;

% Robin BC at x=1
A(end,end) = A(end,end) + alpha;
F(end) = F(end) + alpha*gB;

%% Solve FEM
u_fem = A\F;

%% Analytical solution
syms C
C1 = (alpha*gB - (f/a)*alpha*(1 - exp(-sqrt(a))) - (f/sqrt(a))*exp(-sqrt(a))) ...
      / (alpha*sinh(sqrt(a)) + sqrt(a)*cosh(sqrt(a)));
u_exact = @(x) C1*sinh(sqrt(a)*x) + (f/a)*(1 - exp(-sqrt(a)*x));
u_an = u_exact(x);

%% Plot
plot(x,u_fem,'o-','LineWidth',2); hold on;
plot(x,u_an,'r--','LineWidth',2);
legend('FEM','Analytical');
xlabel('x'); ylabel('u(x)');
title('FEM vs Analytical solution');
grid on;