nn=7; %How many nodes per side

n=nn-2; %interior nodes

h=1/(n+1); %distance 

X=diag(ones(n*n-1,1),1)+diag(ones(n*n-1,1),-1)+diag(ones(n*n-n,1),n)+diag(ones(n*n-n,1),-n)-4*diag(ones(n*n,1),0);
X=-X/(h*h);

for i=1:(n-1)

    X(i*n, i*n+1)=0;X(i*n+1, i*n)=0;
    
end

f=10*ones(n*n,1);

c=X\f;


c2d=zeros(n+2);

for i=1:n
  
      c2d(:,i+1)=[0;c(((i-1)*n+1):(i*n));0];
      
end



[XX,YY]=meshgrid(linspace(0,1,n+2), linspace(0,1,n+2));

         surf(XX,YY,c2d)
            axis([0 1 0 1 -1 1])
