function [phi,exact] = source(X,Y,phi,xs,ys,xs2,ys2,h) 
%% set the point source in a radius of size 2h
%% also return the exact solution 
exact = min(sqrt((X-xs).^2+(Y-ys).^2 ),sqrt((X-xs2).^2+(Y-ys2).^2 )) ;
I = find(exact < 0.1) ;  %% in a fixed radius 
phi(I) = exact(I) ;