% Solve an inhomogeneous ODE by Fourier transform.
function [x,t] = solve_inhomogeneous(FDfunc, f, dt)
% FDfunc = pointer to a function that returns the FFT of some differential
%     operator, D, as a function of an input array of angular frequencies.
%     The ODE is taken to be D x(t) = f(t).
% f = array of N values of the driver f(t).
% x = array of x values given the driver f.
% t = array of times, calculated here for convenient plotting by the end user.

N = numel(f); % find size of time domain
t = ( 0:(N-1) ) * dt; % construct array of times

omega = omega_array(N, dt); % angular frequency (radians per unit time)
fft_D = FDfunc(omega); % FFT of differential operator for all omega values
F = fftn(f); % Fourier transform of the driver
X = F./fft_D; % Fourier transform of the solution
x = ifftn(X); % solution

% I could have squeezed the last 5 lines into one:
% x = ifftn(fftn(f)./FDfunc(omega_array(N, dt)));
% ...but I prefer readable code!
