vvmodeuler2Dneumann(varargin)

function [w1, w2] = vvmodeuler2Dneumann(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 (vector 1by2) - upper bound of spatial (x) variable
%       (Default L = [1,1])
%   T (scalar) - upper bound of time (t) variable
%       (Default T = 10)
%   alpha (scalar) - square root of coefficient of u_{xx} term
%       (Default alpha = 5)
%   m (vector 1by2) - number of discrete spatial intervals
%       (Default m = [100,100])
%   n (scalar) - number of discrete time intervals
%       (Default n = 100)
%
%   w (m+1 x 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{:});
L_x = L(1);
L_y = L(2);
m_x = m(1);
m_y = m(2);

% initialize h, k, lambda, and w
h_x = L_x/m_x;
h_y = L_y/m_y;
k = T./(7500*n);
lambda_x = (alpha.^2).*k./(h_x.^2);
lambda_y = (alpha.^2).*k./(h_y.^2);
w1 = zeros(2*m_y+1,2*m_x+1,n);
w2 = zeros(2*m_y+1,2*m_x+1,n);

% initialize w1 (to sin(pix1/l)*sin(pix2/l)
for i=1:2*m_x+1
    for j=1:2*m_y+1
        w1(j,i,1) = cos(pi.*h_x.*(i-1)./L_x).*cos(pi.*h_y.*(j-1)./L_y);
    end
end

% initialize w2 (to 4*sin(pix1/l)*sin(pix2/l)
for i=1:2*m_x+1
    for j=1:2*m_y+1
        w2(j,i,1) = cos(pi.*h_x.*(i-1)./L_x).*cos(pi.*h_y.*(j-1)./L_y);
    end
end

for t=1:(n-1)
    f1 = laplacian(w1(:,:,t),h_y,h_x);
    y1 = w1(:,:,t) + (k/2)*f1;

    f2 = laplacian(w2(:,:,t),h_y,h_x);
    y2 = w2(:,:,t) + (k/2)*f2;

    % went out to length 2m
    X1 = [sqrt(1./(2*m_y)).*real(fft(y1, 2*m_y,1));zeros(1,201)];

    X2 = [sqrt(1./(2*m_y)).*real(fft(y2, 2*m_y,1));zeros(1,201)];

    Y1 = [sqrt(1./(2*m_x)).*real(fft(X1, 2*m_x,2)) zeros(201,1)];

    Y2 = [sqrt(1./(2*m_x)).*real(fft(X2, 2*m_x,2)) zeros(201,1)];

    % angles divided by m
    a_x = repmat((pi./m_x).*(0:2*m_x),[2*m_y+1,1]);  
    a_y = repmat((pi./m_y).*(0:2*m_y).',[1,2*m_x+1]);   

    % changed c to accomodate for 2-D
    c1 = 1 + lambda_x*lambda_y - 2*(lambda_x.^2)*(lambda_y.^2) + ...
        4*((lambda_x.^2)*(lambda_y.^2).*(cos(a_x).*cos(a_y))) - ...
        ((lambda_x)*(lambda_y).*(cos(a_x).*cos(a_y))) - ...
        2*((lambda_x.^2).*cos(a_x)).*((lambda_y.^2).*cos(a_y));
   
    c2 = 1 + lambda_x*lambda_y - 2*(lambda_x.^2)*(lambda_y.^2) + ...
        4*((lambda_x.^2)*(lambda_y.^2).*(cos(a_x).*cos(a_y))) - ...
        ((lambda_x)*(lambda_y).*(cos(a_x).*cos(a_y))) - ...
        2*((lambda_x.^2).*cos(a_x)).*((lambda_y.^2).*cos(a_y));


    % went out to length 2m
    Z1 = sqrt(1./(2*m_y)).*real(fft(Y1./c1, 2*m_y,1));
    Z1 = Z1([1:2*m_y 1],:);

    Z2 = sqrt(1./(2*m_y)).*real(fft(Y2./c2, 2*m_y,1));
    Z2 = Z2([1:2*m_y 1],:);

    R1 = sqrt(1./(2*m_x)).*real(fft(Z1, 2*m_x,2));
    R1 = R1(:,[1:2*m_x 1]);
    
    R2 = sqrt(1./(2*m_x)).*real(fft(Z2, 2*m_x,2));
    R2 = R2(:,[1:2*m_x 1]);
    
    w1(:,:,t+1) = R1;
    w2(:,:,t+1) = R2;
end

%crop w1 and w2
w1 = w1(1:m_y+1,1:m_x+1,:);
w2 = w2(1:m_y+1,1:m_x+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,100];
end

% check m has valid value
if (numel(m)~=2) || any(m<1) || any(~isequal(round(m),m))
    error('m must be a vector of positive integers.');
end

% get/set alpha
if nargs<3
    alpha = [];
else
    alpha = varargin{3};
end
if isempty(alpha)
    alpha = 5;
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 = 10;
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,1];
end

% check L has valid value
if (numel(L)~=2) || any(L<=0)
    error('L must be a positive vector with 2 elements.');
end

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

%another subfunction
function y = laplacian(x,h_y,h_x)
[m,n] = size(x);
y = zeros(m,n);

y(1,1) = (2*x(1,2)-2*x(1,1))/(h_x.^2)+(2*x(2,1)-2*x(1,1))/(h_y.^2);
for i = 2:n-1
    y(1,i) = (x(1,i+1)-2*x(1,i)+x(1,i-1))/(h_x.^2)+(2*x(2,i)-2*x(1,i))/(h_y.^2);
end
y(1,n) = (2*x(1,n-1)-2*x(1,n))/(h_x.^2)+(2*x(2,n)-2*x(1,n))/(h_y.^2);

for j = 2:m-1
   y(j,1) = (x(j-1,1)+x(j+1,1)-2*x(j,1))/(h_y.^2)+(2*x(j,2)-2*x(j,1))/(h_x.^2);
end

for i = 2:n-1
    for j = 2:m-1
        
    y(j,i) = (x(j-1,i)+x(j+1,i)-2*x(j,i))/(h_y.^2)+(x(j,i-1)+x(j,i+1)-2*x(j,i))/(h_x.^2);
    
    end
end

for j = 2:m-1
   y(j,n) = (x(j-1,n)+x(j+1,n)-2*x(j,n))/(h_y.^2)+(2*x(j,n-1)-2*x(j,n))/(h_x.^2);
end

y(m,1) = (2*x(m,2)-2*x(m,1))/(h_x.^2)+(2*x(m-1,1)-2*x(m,1))/(h_y.^2);
for i = 2:n-1
    y(m,i) = (x(m,i+1)-2*x(m,i)+x(m,i-1))/(h_x.^2)+(2*x(m-1,i)-2*x(m,i))/(h_y.^2);
end
y(m,n) = (2*x(m,n-1)-2*x(m,n))/(h_x.^2)+(2*x(m-1,n)-2*x(m,n))/(h_y.^2);