clear all
% equazione 
%-u''(x)+u(x) = f(x)
%u(a) = u(b)
mrif = 128;
a = -pi;
b = pi;
x = linspace(a,b,mrif+1)';
x = x(1:mrif);
f = @(x) 1./(sin(x)+2);
fhat = fftshift(fft(f(x)))*sqrt((b-a))/mrif;
lambda = i*2*pi*[-mrif/2:mrif/2-1]'/(b-a);
uhatrif = fhat./(1-lambda.^2);
urif = ifft(ifftshift(uhatrif))*mrif/sqrt(b-a);
mrange = [8:8:64];
counter = 0;
for m = mrange
  counter = counter+1;
  x = linspace(a,b,m+1)';
  x = x(1:m);
  fhat = fftshift(fft(f(x)))*sqrt((b-a))/m;
  lambda = i*2*pi*[-m/2:m/2-1]'/(b-a);
  uhat = fhat./(1-lambda.^2);
  u = ifft(ifftshift(uhat))*m/sqrt(b-a);
  errore(counter) = norm(uhatrif-[zeros((mrif-m)/2,1);uhat;zeros((mrif-m)/2,1)]);
end
loglog(mrange,errore,'*',mrange,errore(1)*(mrange/mrange(1)).^(-1),...
mrange,errore(1)*(mrange/mrange(1)).^(-2),...
mrange,errore(1)*(mrange/mrange(1)).^(-4),...
mrange,errore(1)*(mrange/mrange(1)).^(-8))
xlabel('m')
ylabel('errore')
legend('errore','ordine 1','ordine 2','ordine 4','ordine 8','location','NorthWest')