% Poisson Equation via K2D9 for delta(x-1/2,y-1/2), u(unit square)=0, u(m/8,m/8) m=1:7, as functions of n
% David Bloore, August 2011

clc;
clear all;

n=input('Initial n?');
ni=n;

m=input('Number of points along the n axis?');

X1=sparse(zeros(m,7));
Y1=sparse(zeros(m,7));
Z1=sparse(zeros(m,7));

counter=1;

while counter<m+1;
    deltanode=(n^2+1)/2;
    h=1/(n+1);
    I=sparse(eye(n));
    K=sparse(2*I-diag(ones(1,n-1),1)-diag(ones(1,n-1),-1));
    K2D9=kron(K,(-K+6*I))+kron(I,(6*K));
    F=sparse(zeros(n^2,1));
    F(deltanode)=h^2;
    u=K2D9\F;
    u=full(u);
    X1(counter,1:7)=[1/8 1/4 3/8 1/2 5/8 3/4 7/8];
    Y1(counter,1:7)=[n n n n n n n];
    Z1(counter,1:7)=[u((n+1)/8+((n+1)/8-1)*n) u((n+1)/4+((n+1)/4-1)*n) u(3*(n+1)/8+(3*(n+1)/8-1)*n) u(deltanode) u(5*(n+1)/8+(5*(n+1)/8-1)*n) u(n^2-((n+1)/4-1)-((n+1)/4-1)*n) u(7*(n+1)/8+(7*(n+1)/8-1)*n)];
    n=n+8;
    counter=counter+1;
end

plot=plot3(X1,Y1,Z1);
xlim([0 1]);
ylim([ni n-8]);
xlabel('1/8,1/8 1/4,1/4 3/8,3/8 1/2,1/2 5/8,5/8 3/4,3/4 7/8,7/8');
ylabel([sprintf('%d',ni) ' \leq n \leq ' sprintf('%d',n)]);
zlabel('u(1/4,1/4) u(1/2,1/2) u(3/4,3/4)');
title('Poisson Equation via K2D9 for delta(x-1/2,y-1/2), u(unit square)=0, u(m/8,m/8) m=1:7 as functions of n');
set(gca,'XTick',0:1/8:1);
set(gca,'YDir','reverse');
view(-45,20);
% view(-45,25);
grid on;