diff --git a/Sheet4/Ex1.m b/Sheet4/Ex1.m
new file mode 100644
index 0000000..7745237
--- /dev/null
+++ b/Sheet4/Ex1.m
@@ -0,0 +1,50 @@
+%% 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;