clear all
close all
m = 9; % numero dispari, m-1 multiplo di 4
x = linspace(0,1,m)';
epsilon = 1/4;
gepsilon = @(x) (1/2-epsilon <= x & x < 1/2+epsilon) * (-1/(2*epsilon));
h = diff(x);
A = spdiags([[-1./h;0],[1/h(1);1./h(1:m-2)+1./h(2:m-1);1/h(m-1)],[0;-1./h]],[-1,0,1],m,m);
g = zeros(size(x));
n1 = fix((m-1)/4)+1;
n2 = fix((m-1)/4*3)+1;
g(n1) = -1/(2*epsilon) * h(n1)/2;
g(n1+1:n2-1) = -1/(2*epsilon) * (h(n1:n2-2)+h(n1+1:n2-1))/2;
g(n2) = -1/(2*epsilon) * h(n2-1)/2;
% condizioni al bordo (penalizzazione)
M = 2^60;
A(1,1) = M;
g(1) = M*0;
A(m,m) = M;
g(m) = M*0;
uhat = A\g;
xx = linspace(0,1,1001);
idx1 = 0 <= xx & xx < 1/2-epsilon;
idx2 = 1/2-epsilon <= xx & xx < 1/2+epsilon;
idx3 = 1/2+epsilon <= xx & xx <= 1;
uesatta = zeros(size(xx));
uesatta(idx1) = -1/2*xx(idx1);
uesatta(idx2) = 1/(4*epsilon)*(xx(idx2)-1/2).^2+(epsilon-1)/4;
uesatta(idx3) = -1/2*(1-xx(idx3));
plot(x,uhat,'-*',xx,uesatta)
