%We define the boundary gradient of a known harmonic function
  %Construct samples of the input harmonic function along the target grid
    
        %If f is analytic, then u=real(f) is harmonic, and
        %grad u=(real(f'),-imag(f')), which simplifies computations.
        n=5
         display('Try n = 6! ')
        hh(:,:)=real((x(:,:)+i*y(:,:)).^n);
        grdh(:,1)=real(n*(x2(1:N)+i*y2(1:N)).^(n-1));
        grdh(:,2)=-imag(n*(x2(1:N)+i*y2(1:N)).^(n-1));
    hnu=(sum((nrm_vec.*grdh)'))';
    ww=Dm1\(hnu);
    %Check that hnu is in the range of the normal derivative of the single
    %layer
    display('The projection of the data onto the nullspace of the adjoint is: ')
    u(:,N)'*hnu
    %Check to see that we are getting the solution perpindicular to the
    %null-space.
    display('The projection of the density onto the nullspace of the double layer is: ')
    v(:,N)'*ww
    ln=length(ww);
    figure, plot(1:ln,ww,'r',1:ln,hnu,'b')
    title('The density on the boundary (red) and the normal derivative (blue)')
     for j=1:Nr
        for k=1:Nth
            hhc(j,k)=0;
            for l=1:N
            hhc(j,k)=sl(j,k,l)*ww(l)+hhc(j,k);
            end
        end
     end
     figure, surf(x,y,hhc-hh)
     shading interp
     title('The error')
     ax=gca;
     
     xl=ax.XLim;
     yl=ax.YLim;
     axlim=max(xl(2),yl(2));
     ax.XLim=[-axlim,axlim];
     ax.YLim=[-axlim,axlim];
     figure, surf(x,y,hhc)
     title('The computed solution')
     shading interp
     ax=gca;
     
     xl=ax.XLim;
     yl=ax.YLim;
     axlim=max(xl(2),yl(2));
     ax.XLim=[-axlim,axlim];
     ax.YLim=[-axlim,axlim];
    
     %
     %
     display('The maximum error is:')
     max(max(abs(hh-hhc)))
    