%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% second spatial derivative (No flux BCs)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [result] = diff2d_conic(xv,Nx,Ny,dyy,dxx,dyyy)


%sides
i=1; j=2:Nx-1;
result(i,j)=dyy(j).*(xv(i+1,j)+xv(i+1,j)-2*xv(i,j))+ ...
            dxx*(xv(i,j+1)+xv(i,j-1)-2*xv(i,j));
result(Ny,j)=result(1,j);
i=2:Ny-1; j=1;
result(i,j)=dyy(j).*(xv(i+1,j)+xv(i-1,j)-2*xv(i,j))+ ...
            dxx*(xv(i,j+1)+xv(i,j+1)-2*xv(i,j));
i=2:Ny-1; j=Nx;
result(i,j)=dyy(j).*(xv(i+1,j)+xv(i-1,j)-2*xv(i,j))+ ...
            dxx*(xv(i,j-1)+xv(i,j-1)-2*xv(i,j));

%corners
i=1; j=1;
result(i,j)=dyy(j).*(xv(i+1,j)+xv(i+1,j)-2*xv(i,j))+ ...
            dxx*(xv(i,j+1)+xv(i,j+1)-2*xv(i,j));
i=1; j=Nx;
result(i,j)=dyy(j).*(xv(i+1,j)+xv(i+1,j)-2*xv(i,j))+ ...
            dxx*(xv(i,j-1)+xv(i,j-1)-2*xv(i,j));
result(Ny,1)= result(1,1);
result(Ny,Nx)= result(1,Nx);

%middle
i=2:Ny-1; j=2:Nx-1;
result(i,j)=dyyy.*(xv(i+1,j)+xv(i-1,j)-2*xv(i,j))+ ...
            dxx*(xv(i,j+1)+xv(i,j-1)-2*xv(i,j));


return

