cnfftneumann(varargin)

function w = cnfftneumann(varargin)
% crankNicolson: uses Crank-Nicolson algorithm to approximate the solution
%   to the parabolic PDE:
%       u_{t}(x,t) - alpha^2 u_{xx}(x,t) = 0, 0<x<L, 0<t<T
%   subject to the boundary conditions
%       u(0,t) = u(L,t) = 0, 0<t<T
%   and the initial conditions
%       u(x,0) = f(x), 0<=x<=L
%
% arguments:
%   L (scalar) - upper bound of spatial (x) variable
%       (Default L = 1)
%   T (scalar) - upper bound of time (t) variable
%       (Default T = 1)
%   alpha (scalar) - square root of coefficient of u_{xx} term
%       (Default alpha = 1)
%   m (scalar) - number of discrete spatial intervals
%       (Default m = 100)
%   n (scalar) - number of discrete time intervals
%       (Default n = 100)
%
%   w (m+1 x n) - approximation to u(x,t) at discrete space/time positions
%

% author: Troy J. Winkstern
% email: [e-mail]
% date: 30 Jan 2011

% parse input arguments
[L,T,alpha,m,n] = parseInputs(varargin{:});

% initialize h, k, lambda, and w
h = L./m;
k = T./n;
lambda = (alpha.^2).*k./(h.^2);
%w = zeros(m+1,n);
w = zeros(2*m+1,n);

% initialize first column of w
w(:,1) = 5*cos(pi.*h.*(0:2*m).');

for i=1:(n-1)
   y = zeros(2*m+1,1);
   
   y(1) = lambda*w(2,i) + (1-lambda)*w(1,i);
   for j = 2:(2*m)
       y(j) = (lambda/2)*w(j-1,i) + (1-lambda)*w(j,i) + (lambda/2)*w(j+1,i);
   end
   y(2*m+1) = lambda*w(2*m,i) + (1-lambda)*w(2*m+1,i);
   
   % went out to length 2m
   X = [sqrt(1./(2*m)).*real(fft(y,2*m));0];
   
   % angles divided by m
   a = (pi./m).*(0:2*m).';
   
   % changed lambda to negative
   c = -lambda*(cos(a)) + 1 + lambda;

   % went out to length 2m
   Z = sqrt(1./(2*m)).*real(fft(X./c, 2*m));
   w(1:2*m,i+1) = Z(1:2*m);
   w(2*m+1,i+1) = w(1,i+1);
end

% crop w
w = w(1:m+1,:);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction parseInputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [L,T,alpha,m,n] = parseInputs(varargin)

% check number of input arguments
nargs = length(varargin);
error(nargchk(0,5,nargs));

% get/set n
if nargs<5
    n = [];
else
    n = varargin{5};
end
if isempty(n)
    n = 100;
end

% check n has valid value
if (numel(n)>1) || (n<1) || ~isequal(round(n),n)
    error('n must be a positive integer.');
end

% get/set m
if nargs<4
    m = [];
else
    m = varargin{4};
end
if isempty(m)
    m = 100;
end

% check m has valid value
if (numel(m)>1) || (m<1) || ~isequal(round(m),m)
    error('m must be a positive integer.');
end

% get/set alpha
if nargs<3
    alpha = [];
else
    alpha = varargin{3};
end
if isempty(alpha)
    alpha = 1;
end

% check alpha has valid value
if numel(alpha)>1
    error('alpha must be a real scalar.');
end

% get/set T
if nargs<2
    T = [];
else
    T = varargin{2};
end
if isempty(T)
    T = 1;
end

% check T has valid value
if (numel(T)>1) || (T<=0)
    error('T must be a positive real scalar.');
end

% get/set L
if nargs<1
    L = [];
else
    L = varargin{1};
end
if isempty(L)
    L = 1;
end

% check L has valid value
if (numel(L)>1) || (L<=0)
    error('L must be a positive real scalar.');
end

% warning if L not integer
if ~isequal(round(L),L)
    warning('L should be an integer so that the inital condition takes on the value of 0 for t = 0 and t = T.');
end