%
% equazione 
% -u''(x)+u(x) = g(x), -pi <= x < pi
% u(a) = u(b)
%
mrif = 256;
m = mrif;
a = -pi;
b = pi;
x = linspace (a, b, m + 1)';
x = x(1:m);
g = @(x) 1 ./ (sin (x) + 2);
ghat = fftshift (fft (g(x))) * sqrt (b - a) / m;
lambda = 1i * 2 * pi * (-m/2:m/2-1)' / (b - a);
uhat = ghat./(1-lambda.^2);
u = ifft (ifftshift (uhat)) * m / sqrt(b - a);
uhatrif = uhat;
urif = u;
mrange = [8:8:mrif-16];
counter = 0;
for m = mrange
  counter = counter+1;
  x = linspace (a, b, m + 1)';
  x = x(1:m);
  ghat = fftshift (fft (g(x))) * sqrt (b - a) / m;
  lambda = 1i * 2 * pi * (-m/2:m/2-1)' / (b - a);
  uhat = ghat ./ (1 - lambda .^ 2);
  u = ifft (ifftshift (uhat)) * m / sqrt (b - a);
  uhat = [zeros((mrif - m) / 2, 1); uhat; zeros((mrif - m) / 2, 1)];
% Parseval: norma due della differenza dei coefficienti  
  errore(counter) = norm (uhatrif-uhat);
end
loglog(mrange, errore(1) * (mrange / mrange(1)) .^ (-1),'k',...
mrange,errore(1) * (mrange / mrange(1)) .^ (-2),'r',...
mrange,errore(1) * (mrange / mrange(1)) .^ (-4),'m',...
mrange,errore(1) * (mrange / mrange(1)) .^ (-8),'g',mrange, errore, '*b')
xlabel('m')
ylabel('errore in norma-2')
legend('ordine 1','ordine 2','ordine 4','ordine 8','errore',...
'location','NorthWest')
