function A_1_4(n)
% Anzahl der Punkte als input.

h=2/(n-1); % Schrittweite - beachte Intervallbreite von 2
x = (-1+h:h:1-h)'; % äquidistantes Gitter der inneren Punkte
x_zwischen = -1+h/2:h:1-h/2; % Zwischenpunkte
f=ones(size(x)); % rechte Seite
% definiere function handle für k
k = @(x) 1+(x>0);
% evaluiere k an den Zwischenpunkten
k_eval=k(x_zwischen);
k_links = k_eval(1:end-1)';
k_rechts = k_eval(2:end)';  
% baue Matrix A
A = sparse(diag(-k_links(2:end),-1) + diag(-k_rechts(1:end-1),1) + diag(k_links+k_rechts,0));
% alternative über spdiags
%Atest = spdiags([-[k_links(2:end);0], [k_links+k_rechts], -[0;k_rechts(1:end-1)]],-1:1,length(x),length(x));
% löse LGS
u=(1/h^2*A)\f;
% Ränder
u = [0;u;0];
% x inklusive Randpunkte
x_full = [-1;x;1];
% function handle der schwachen Lösung - hier mit u_weak(0) = 1/3 - kann
% man auch weglassen, dann wird die schwache Lösung in der 0 auf 0 gesetzt.
u_weak = @(x) (-x.^2./2 - x./6 + 1/3).*(x <= 0) + (-x.^2./4 - x./12 + 1/3).*(x > 0);
% plot
figure; 
plot(x_full,u,'k-');
% weak solution
hold on
plot(x_full,u_weak(x_full), 'r-');
title('FD Approximation')
xlabel('x')
ylabel('u(x)','interpreter','latex')
legend('FD','weak solution')

% nicht C^2, wird deutlich mit:
%figure; plot(diff(u));