%	risoluzione mediante la trasformata di Fourier della ODE
%		y''+6y'+5y = f(x)
%	per alcune f(x)

% 	Giuseppe Rodriguez,  27-09-05
% 	Email: rodriguez@unica.it

A = 40;
N = 256;
h = A/N;
W = N/A;

x = [0:N-1]'/N*A-A/2;
w = [0:N-1]'/N*W-W/2;

% esempio 1
g = exp(-x.*x);
f = exp(-x.*x) .* (4*x.*x-12*x+3);

% esempio 2
%g = sin(x)./x;
%g(x==0) = 1;
%f = 2*((1-3*x).*(sin(x)-x.*cos(x))./x./x./x + 2*sin(x)./x);
%f(x==0) = 14/3;

F = A/N * fft(fftshift(f));
F = fftshift(F);

Y = F ./ (5+12*pi*i*w-4*pi*pi*w.*w);
y = 1/A*N * ifft(fftshift(Y));
y = fftshift(y);
norm(imag(y))
y = real(y);

figure(1)
subplot(2,1,1)
plot(x,real(F));
title('funzione F')
xlabel('parte reale')
subplot(2,1,2)
plot(x,imag(F));
xlabel('parte immaginaria')
figure(2)
plot(x,y,x,g,'--')
title('soluzione')
figure(3)
plot(x,y-g)
title('errore')

