clear all
close all
m = 80;
epsilon = 1/4;
gepsilon = @(x) (1/2-epsilon <= x & x < 1/2+epsilon) * (-1/(2*epsilon));
xl = linspace (0, 1/2-epsilon, round(m/3))';
x = xl;
n1 = length (x);
xm = linspace (1/2-epsilon, 1/2+epsilon, round(m/3))';
x = [x;xm(2:end)];
n2 = length (x);
xr = linspace (1/2+epsilon, 1, round(m/3))';
x = [x;xr(2:end)];
m = length (x);
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));
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;
idx1 = 0 <= x & x < 1/2-epsilon;
idx2 = 1/2-epsilon <= x & x < 1/2+epsilon;
idx3 = 1/2+epsilon <= x & x <= 1;
uesatta = zeros(size(x));
uesatta(idx1) = -1/2*x(idx1);
uesatta(idx2) = 1/(4*epsilon)*(x(idx2)-1/2).^2+(epsilon-1)/4;
uesatta(idx3) = -1/2*(1-x(idx3));
plot(x,uhat,'*',x,uesatta)
norm(uesatta-uhat,inf)
mm = 1001;
xx = linspace(0,1,mm)';
uuhat = interp1(x,uhat,xx,'linear');
idx1 = 0 <= xx & xx < 1/2-epsilon;
idx2 = 1/2-epsilon <= xx & xx < 1/2+epsilon;
idx3 = 1/2+epsilon <= xx & xx <= 1;
uuesatta = zeros(size(xx));
uuesatta(idx1) = -1/2*xx(idx1);
uuesatta(idx2) = 1/(4*epsilon)*(xx(idx2)-1/2).^2+(epsilon-1)/4;
uuesatta(idx3) = -1/2*(1-xx(idx3));
norm(uuhat-uuesatta)/sqrt(mm-1)
