function  dphidr= phiodek(r,phi,ri,k,v,F1,F2)

x1=interp1(ri,F1,r,'spline');
x2=interp1(ri,F2,r,'spline');

dphidr = [ 
           phi(3) 
           phi(4)
           -phi(3)./r+((2+k^2*(1-2*v)-2*v)/2/(1-v)).*phi(1)./r.^2-(k*(3-4*v)/2/(1-v)).*phi(2)./r.^2+k/2/(1-v).*phi(4)./r+x1
           -phi(4)./r-(k*(3-4*v)/(1-2*v)).*phi(1)./r.^2+((1+2*k^2*(1-v)-2*v)/(1-2*v)).*phi(2)./r.^2-k/(1-2*v).*phi(3)./r+x2
                   ];