There are a few issues with your code. The first is the use of linspace. It includes both endpoints of the interval, thus both $0$ and $4\pi$ appear. This is not what you want to happen with the discretization for the purpose of Fourier transform. E.g., linspace(0,2*pi,4) returns the vector $0,2\pi/3, 4\pi/3, 2\pi$ but what we actually want is $0,\pi/2, \pi , 3\pi/2$. Solution: introduce the step $dx=2\pi/N$ and create the vector a+[0:N-1]*dx.
Second, the correct version of $2\pi i \xi$ in the discrete setting is not obvious, due to multiple ways to enumerate the terms of Fourier series. In the code below this role is played by vector $k$. I adapted it from Finding Derivatives using Fourier Spectral Methods.
a = 0;
b = 4*pi;
N = 4096;
dx = (b-a)/N;
t = a + dx*(0:N-1);
f = sin(t);
fftx = fft(f);
k = 2*pi/(b-a)*[0:N/2-1, 0, -N/2+1:-1];
dffft = 1i*k.*fftx;
df = ifft(dffft);
plot(t,df)
Why the strange $k$ vector? Part of it is scaling: the "wider" the function, the more narrow is its Fourier transform (it scales in the opposite direction). This is why we divide by $(b-a)$ there. To understand the weird [0:N/2-1, 0, -N/2+1:-1], think of what the components of fft are: they are the coefficients of the complex polynomial that matches $f$. For a simple example, take t=0:pi/4:(2*pi-pi/4); that is, $N=8$ here. Applying fft to
f = 3+7*exp(1i*t)+8*exp(3i*t)+9*exp(5i*t)+10*exp(7i*t)
we get 8*[3 7 0 8 0 9 0 10] which is, as expected, the coefficients of exponentials in $f$ (multiplied by $8$, which isn't important). Multiplying the coefficients by 1i*[0:7] gives us the Fourier coefficients of $f'$, and ifft returns $f'$ itself. Sweet. As long as we work with polynomials of $e^{it}$, all this works fine.
Trouble arises when our $f$ has negative powers of $e^{it}$, for example $\sin t = \frac{1}{2}(e^{it}-e^{-it})$. Where does $e^{-it}$ go under the fft? Same place where $e^{7it}$ goes, since $e^{-it} \equiv e^{7it}$ on our grid. But... then our algorithm will multiply it by $7i$, not by $-i$! Most assuredly,
$$(\sin t)' \ne \frac{1}{2}(ie^{it}-7ie^{-it})$$
This is the problem that [0:N/2-1, 0, -N/2+1:-1] solves. It simply recognizes that $e^{7it}$ should be treated as $e^{-it}$ for the purpose of taking the derivative: that is, it should be multiplied by $-i$, not by $7i$.
With $N=8$, the multiplying vector is [0 1 2 3 0 -3 -2 -1]. It
- kills the constant term,
- multiplies $e^{ikt}$ component by $k$ for $k=1,2,3$,
- kills the $e^{4it}$ component (What else to do with it? We don't know if it's $e^{4it}$ or $e^{-4it}$. Note that $\sin 4t $ will have zero DFT at this scale).
- multiplies $e^{ikt} = e^{i(k-8)t}$ component by $k-8$ for $k=5,6,7$.
It usually helps to follow Fourier's prescription rather than jump to the solution. For Laplace's equation, the separated solutions $X(x)Y(y)$ satisfy
$$
X''Y+XY'' = 0, \\
\frac{X''}{X} = -\frac{Y''}{Y} = \lambda,
$$
where $\lambda$ is a separation parameter. So the Fourier ODEs are
$$
\begin{align}
X''=\lambda X, &\;\;\;\;\; Y''=-\lambda Y \\
0 < x < a, & \;\;\;\;\;0 < y < \infty \\
X(0) = 0, & \;\;\;\;\; Y'(0) = 0.
\end{align}
$$
Fourier always assumed the separated solutions were bounded on any unbounded domain. The conditions for $X$ do not determine parameters $\lambda$; however the boundedness of $Y$ in conjunction with $Y'(0)=0$ forces $\lambda \ge 0$. (If $\lambda < 0$ then $Y(y)=\cosh(\sqrt{\lambda}y)$ is unbounded.) Therefore, $\lambda \ge 0$ is assumed, and the separated solutions for $\lambda=\mu^{2}$ satisfying the above conditions are
$$
X_{\mu}(x)Y_{\mu}(y) = C(\mu)\sinh(\mu x)\cos(\mu y).
$$
A sum of these must be an integral sum:
$$
u(x,y) = \int_{0}^{\infty}C(\mu)\sinh(\mu x)\cos(\mu y)dy.
$$
The coefficient function $C(\mu)$ is the only missing component, and must be chosen so that
$$
g(y) = u(a,y) = \int_{0}^{\infty}C(\mu)\sinh(\mu a)\cos(\mu y)dy.
$$
Writing $g$ as a Fourier cosine transform gives
$$
g(y) = \frac{2}{\pi}\int_{0}^{\infty}\left(\int_{0}^{\infty}\cos(\mu u)g(u)du\right)\cos(\mu y)d\mu,
$$
which leads to
$$
C(\mu)\sinh(\mu a) = \frac{2}{\pi}\int_{0}^{\infty}\cos(\mu u)g(u)du \\
C(\mu) = \frac{2}{\pi\sinh(\mu a)}\int_{0}^{\infty}\cos(\mu u)g(u)du.
$$
The final Fourier solution is
$$
u(x,y) = \frac{2}{\pi}\int_{0}^{\infty}\left(\frac{1}{\sinh(\mu a)}
\int_{0}^{\infty}\cos(\mu u)g(u)du\right)\sinh(\mu x) \cos(\mu y)d\mu.
$$
Sanity Check: Note that the expression in parentheses is a function of $\mu$ only--it is $C(\mu)$. Visual inspection of the final, proposed solution shows that $u(a,y)$ is the inverse Fourier cosine transform of the Fourier cosine transform of $g$. So $u(a,y)=g(y)$ under suitable smoothness conditions on $g$. Cleary $u(0,y)=0$. Under suitable assumptions to ensure the convergence of the differentiated integral expression, $u_{y}(x,0)=0$ follows. All of the required conditions check, and this does give a solution of Laplace's equation given sufficient convergence of the different derived series.
Best Answer
The Neumann boundary condition asks that we solve a differential equation subject to the condition that the derivative of the function involved takes on some specific value at the boundary of the domain over which we are solving. It follows that we require the derivative to exist at the boundary but other than that it does not effect our ability to represent the derivative as a cosine transform.
Assuming that the Fourier transform of the function exists, the following shows when the Fourier transform representation of the derivative reduces to a cosine transform; it does not require the function to be periodic.
If you define: \begin{align}\hat{F}(\xi)&=\int_{-\infty}^{\infty}f(x)e^{-i2\pi x \xi}{dx},\tag{1}\\ f(x)&=\int_{-\infty}^{\infty}\hat{F}(\xi)e^{i2\pi x \xi}{d\xi},\end{align} then: \begin{aligned} \dfrac{d}{dx}f(x)&=\int_{-\infty}^{\infty}\left(i2\pi\xi\right)\hat{F}(\xi)e^{i2\pi x \xi}{d\xi},\\ \dfrac{d}{dx}f(x)&=\int_{-\infty}^{\infty}\left(i2\pi\xi\right)\hat{F}(\xi)\cos(2\pi x \xi){d\xi}-\int_{-\infty}^{\infty}\left(2\pi\xi\right)\hat{F}(\xi)\sin(2\pi x \xi){d\xi},\end{aligned}
so to reduce the Fourier transform representation of the derivative to the cosine transform only we would require: $$\int_{-\infty}^{\infty}\left(2\pi\xi\right)\hat{F}(\xi)\sin(2\pi x \xi){d\xi}=0$$ and as $\sin$ is an odd function this integral vanishes when $\left(2\pi\xi\right)\hat{F}(\xi)$ is an even function which occurs when: $$\xi\hat{F}(\xi)=-\xi\hat{F}(-\xi),$$ $$\hat{F}(\xi)=-\hat{F}(-\xi),\tag{2a}$$ that is when $\hat{F}(\xi),$ is odd, and we see from $(1)$ that this requires $$f(x)=-f(-x)\tag{2b}$$ i.e. it requires $f(x)$ to be odd. But if $(2)$ holds then $(1)$ can be written: \begin{align} \hat{F}(\xi)&=-i\int_{-\infty}^{\infty}f(x)\sin(2\pi x \xi){dx},\tag{3}\\ f(x)&=i\int_{-\infty}^{\infty}\hat{F}(\xi)\sin(2\pi x \xi){d\xi}, \end{align}
So to summarise, assuming $(1)$ exists and $(2)$ holds we may define the derivative by the following cosine transform: $$\dfrac{d}{dx}f(x)=i2\pi\int_{-\infty}^{\infty}\xi\hat{F}(\xi)\cos(2\pi x \xi){d\xi}$$ which we see would also result from applying differentiation under the integral sign to $(3)$. With that said, you have not specified what it is that you want to cosine transform in order to find the derivative; if this equation holds: $$\dfrac{d}{dx}f(x)=\int_{-\infty}^{\infty}\left(i2\pi\xi\right)\hat{F}(\xi)e^{i2\pi x \xi}{d\xi},$$ then we could simply define: $$\hat{F}(\xi)=\hat{G}(\xi)\cos(2\pi x \xi)e^{-i2\pi x \xi},$$ to obtain a cosine transform representation of the derivative but this is being a bit pedantic; some clarification of what it is you hope to cosine transform would help to clear that up.