function fig111

clf

nr=50;
rMax=10;
alpha2=0.1; z1=1; z2=-2; alpha1=-alpha2*z2/z1;
k=sqrt(alpha1*z1^2+alpha2*z2^2);
lambda=0.5*(alpha1*z1^3+alpha2*z2^3);
ep=0.1;

r=linspace(1,rMax,nr);
gamma=lambda*exp(2*k)/(2*k)/(1+k)^2;
a=gamma*( (k-1)*exp(2*k)*expint(3*k) + (k+1)*expint(k) )/(k*(1+k));

for ir=1:nr
	ar=k*r(ir);
	phi0(ir)=exp(k-ar)/r(ir)/(1+k);
        phi1(ir) = a*exp(-ar)/r(ir) + gamma*( exp(ar)*expint(3*ar) - exp(-ar)*expint(ar) )/(ar);
end

phi=ep*(phi0+ep*phi1);

% note rM>rMax
rM=20;

% get(gcf)
set(gcf,'Position', [1653 944 618 232]);

%  calculate solution using a MATLAB routine
sol0 = bvpinit(linspace(1,rM,1000),@mat4init);
sol = bvp4c(@rhs,@bcs,sol0)

%  plot 
hold on
xint = linspace(1,rM);
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:),'Linewidth',1.2)
box on

plot(r,ep*phi0,'--r','Linewidth',1.2)
%plot(r,phi,'k')

box on
grid on
axis([0 rMax 0 0.06])

xlabel('r-axis','FontSize',14,'FontWeight','bold')
ylabel('Potential','FontSize',14,'FontWeight','bold')

set(gca,'FontSize',14);
legend(' \phi',' \epsilon\phi_0','Location','NorthEast');
%legend(' \phi',' \phi_0',' \phi_0 + \epsilon \phi_1','Location','NorthEast');
set(findobj(gcf,'tag','legend'),'FontSize',16); 

hold off


%  define RHS of equations
function dy=rhs(x,y)
dy=zeros(2,1);
dy(1) = y(2);
dy(2) = p(y(1)) - 2*y(2)/x;

%  define BCs
function res = bcs(ya,yb)
ep=0.1;
res = [  ya(2) + ep
         yb(1)  ];

function yinit = mat4init(x)
ep=0.1; 
alpha2=0.1; z1=1; z2=-2; alpha1=-alpha2*z2/z1;
k=sqrt(alpha1*z1^2+alpha2*z2^2);
c = ep*exp(k*(1-x))/((1+k)*x);
yinit = c* [  1
          -k-1/x ];


function g=p(y)
alpha2=0.1; z1=1; z2=-2; alpha1=-alpha2*z2/z1;
g= - alpha1*z1*exp(-z1*y) - alpha2*z2*exp(-z2*y) ;