% Solves Poisson Equation via K2D9 for delta(x-1/2,y-1/2), u(unit square)=0 -- Generates Mesh Plot, frames correspond to increases of N by 2
% David Bloore, August 2011

clc;
clear all;

% input for initial N includes a test to ensure it is an odd integer

boolean=0;

while boolean<1;
    n=input('Initial N?');
    nn=floor(n);
    if nn>1;
        if nn==n;
            nnn=nn^2/2;
            nnnn=floor(nnn);
            if nnnn<nnn
                boolean=1;
            else
                clc;
                display('Please enter an odd integer');
            end
        else
            clc;
            display('Please enter an odd integer');
        end
    else
        clc;
        display('Please enter an odd integer');
    end
end

nf=input('Number of frames?');
nff=2*nf+n+1;

% blank matrices to initate the plot

X=zeros(2);
Y=X;
Z=X'*X;

% the plot

z=mesh(X,Y,Z);
view(-30,65);
title('Solution of Poisson Equation via K2D9 with f=delta(x-1/2,y-1/2)');
xlabel('0 \leq x \leq 1');
ylabel('0 \leq y \leq 1');
zlabel('z');
xlim([0 1]);
ylim([0 1]);
set(gcf,'Position',[0 356 1440 900]);

while n<nff;
    
    % global scalar variables constructed from n
    
    deltanode=(n^2+1)/2;
    h=1/(n+1);
    
    % construction of K
    
    I=sparse(eye(n));
    K=sparse(2*I-diag(ones(1,n-1),1)-diag(ones(1,n-1),-1));
    
    % construction of K2D9
    
    K2D9=kron(K,(-K+6*I))+kron(I,(6*K));
    
    F=sparse(zeros(n^2,1));
    F(deltanode)=h^2;
    
    % solution of system
    
    u=K2D9\F;
    
    % construction of matrices for the plot
    
    Z=sparse(zeros(n));
    
    counter1=1;
    while counter1<n^2+1;
        Z(counter1)=u(counter1);
        counter1=counter1+1;
    end
    
    X=linspace(h,1-h,n);
    Y=linspace(h,1-h,n);
    
    set(z,'XData',X);
    set(z,'YData',Y);
    set(z,'ZData',Z);
    
    drawnow
    n=n+2;
end