Solving 1D Boundary Value Problem using Discrete Fourier Transform. I am trying to solve it in Matlab, where did I mess up?

Question

I have a basic steady-state Boundary Value Problem (BVP) that I am trying to solve discretely by using Fourier Transforms. The setup is such that: $$\cases{u(x_i)_{xx} = f(x_i),f(x_i) = 2, & $x\in [0,1]$; \\u(x_{i=0}) = 0, u(x_{i=N}) = 1}$$ and I worked out the following finite difference formulas set up for the BCs: $$u(x_1)_{xx}=\frac{u(x_{2})+u({x_0})-2u(x_1)}{h^2} \rightarrow u(x_1)_{xx}-\frac{u(x_0)}{h^2} = \frac{u(x_{2})-2u(x_1)}{h^2} \rightarrow u(x_1)_{xx} = 2-\frac{1}{h^2}$$ $$\boxed{f(x_1) = 2-\frac{1}{h^2}}$$ $$u(x_{N-1})_{xx}=\frac{u(x_{N})+u(x_{N-2})-2u(x_{N-1})}{h^2} \rightarrow u(x_{N-1})_{xx}-\frac{u(x_N)}{h^2} = \frac{u(x_{N-2})-2u(x_{N-1})}{h^2} \rightarrow u(x_1)_{xx} = 2-\frac{2}{h^2}$$ $$\boxed{f(x_{N-1}) = 2-\frac{2}{h^2}}$$

Code for solving via Fourier Transform below:

% 1-D Poisson Equation: u_xx = f(x)
% u_xx = 2, 0<x<1
% u(0) = 1, u(1) = 2

L = 1;
Ne = 16; % total number of elements. Is a power of 2
h = 2*L/Ne; % unit spacing
m = Ne/2 -1; % number of interior elements

x = h*(1:m);
f = 2*ones(1,length(x));

% Correct for Dirichlet BCs
u0 = 1; 
u1 = 2; 
f(1) = f(1) - u0/(h^2);%f(1,:) = f(1,:) - u0y/h^2;
f(m) = f(m) - u1/(h^2);%f(m,:) = f(m,:) - u1y/h^2;


%odd extension of f in x
g = [0, f, 0, f(m:-1:1)]; % negative reflection across y-axis i.e. odd extension
ghat = fft(g);
I = 0:Ne-1;
mu = (4/h^2)*(sin(I*pi/Ne).^2); % Fourier Transform(d^2/dx^2 u(x)) 
                                % = (e^{-jw}+e^{jw}-2)uhat(e^jw)/h^2
                                % = (cos(w)-jsin(w) + cos(w)+jsin(w) -2)uhat(e^jw)/h^2
                                % = 2(cos(w)-1)uhat(e^jw)/h^2 = 4sin^2(w/2)uhat(e^jw)/h^2
                                % mu = 4sin^2(w/2)/h^2
mu(1) = 1; % avoid division by 0 
v = real(ifft(-ghat./mu));
u = v(1:(m+2)); % Extract u(i,j) from its odd extension
uex = 1+x.^2;

figure(1)
plot(x,uex);
hold on
plot(x,u(2:(m+1)));
hold off
legend('uex','u')

The solution is supposed to be $u(x) = 1+x^2$ yet I get looks more like 22+3x^2 as we can see in the plot below. So it's too high up and too steep.

I don't know what is missing to make this work out. Is the odd extension wrong or is it something else?

Disclaimer: This is not my code, this was code designed to solve the 2D case, that I then modified for a 1D case. The original code was found here: LINK

EDIT: I have added an additional plot for both g and the Fourier Transform of g.

I think this is supposed to be negative on the right half for the odd extension. Would it be positive if it was an even extension or how does that work? Literally the extension is what I'm not okay on at this point.

Answer 1

I found the issue. The odd extension was supposed to be negative but I had it as positive. The fixed line of code should be: g = [0, f, 0, -f(m:-1:1)];

The full program should be

% 1-D Poisson Equation: u_xx = f(x)
% u_xx = 2, 0<x<1
% u(0) = 1, u(1) = 2

L = 1;
Ne = 16; % total number of elements. Is a power of 2
h = 2*L/Ne; % unit spacing
m = Ne/2 -1; % number of interior elements

x = h*(1:m);
f = 2*ones(1,length(x));

% Correct for Dirichlet BCs
u0 = 1; 
u1 = 2; 
f(1) = f(1) - u0/(h^2);%f(1,:) = f(1,:) - u0y/h^2;
f(m) = f(m) - u1/(h^2);%f(m,:) = f(m,:) - u1y/h^2;


%odd extension of f in x
g = [0, f, 0, -f(m:-1:1)]; % negative reflection across y-axis i.e. odd extension
ghat = fft(g);
I = 0:Ne-1;
mu = (4/h^2)*(sin(I*pi/Ne).^2); % Fourier Transform(d^2/dx^2 u(x)) 
                                % = (e^{-jw}+e^{jw}-2)uhat(e^jw)/h^2
                                % = (cos(w)-jsin(w) + cos(w)+jsin(w) -2)uhat(e^jw)/h^2
                                % = 2(cos(w)-1)uhat(e^jw)/h^2 = 4sin^2(w/2)uhat(e^jw)/h^2
                                % mu = 4sin^2(w/2)/h^2
mu(1) = 1; % avoid division by 0 
v = real(ifft(-ghat./mu));
u = v(1:(m+2)); % Extract u(i,j) from its odd extension
uex = 1+x.^2;

figure(1)
plot(x,uex);
hold on
plot(x,u(2:(m+1)));
hold off
legend('uex','u')
```