function  dphidr= phiode0(r,phi,ri,F)

x=interp1(ri,F,r,'spline');

  dphidr= [ 
           phi(2) 
           -phi(2)./r+phi(1)./r.^2+x                   
            ];