function pile21

clear all;   
close all;

% global  y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 m11 c1 k1 P0 M0

r0=0.18;l=10; Ep=2.7e5;rm=20*r0;nth=100;v= 0.2;M0= 0.*r0;Ap=pi.*(r0.^2);
N=100;nn=-0.5;L=l;Ip=(pi./4).*r0^4;AA=600;P0=30;A=Ep.*Ip;L=l;nn=-0.5;
z=zeros(1,N);tol=0.001;   %tolerance
z0=0; z11=3;z22=4; z33=8 ;z44=l; n1=N; n2=2*N; n3=3*N; n4=4*N;dz=L/(4*N);
eps=zeros(4,N);
phi=zeros(4,N);
r = logspace(log10(r0),log10(rm),N);
z1=linspace(0,z11,N);
z2=linspace(z11,z22,N);
z3=linspace(z22,z33,N);
z4=linspace(z33,z44,N);
rt=logspace(log10(r0),log10(rm),N);
%coefficients for eps 
eps1=0.000001*ones(1,N);
eps2=0.000001*ones(1,N);
eps3=0.000001*ones(1,N);
eps4=0.000001*ones(1,N);

t1=dfdx(eps1(1,:),z1,N);
t2=dfdx(eps2(1,:),z2,N);
t3=dfdx(eps3(1,:),z3,N);
t4=dfdx(eps4(1,:),z4,N);
rt=r;
%phi(1)=phi_r,phi(2)=diff(phi_r);
y1new=ones(1,N); y2new=ones(1,N);....
y1old=0.01*ones(1,N);y2old=0.01*ones(1,N);

while(norm((y1old-y1new)./y1old))>tol||(norm((y2old-y2new)./y2old))>tol....    
y1=y1new;y2=y2new; [1/sqrt(N)*norm((y1old-y1new)./y1old),1/sqrt(N)*...
norm((y2old-y2new)./y2old)];
%--------------------------------------------------------------------------
solinit=bvpinit(logspace(log10(r0),log10(rm),N),[0.0001 1]);
sol = bvp4c(@(r,phi)phiode(r,phi,rt,y1,y2),@phibc,solinit);
phi = deval(sol,r);
%--------------------------------------------------------------------------
%coefficients for eps1 
C1=0;k1=0;m1=0;s1=0;f1=0;z0=0;

    mu1=Gs11(z0,z1,rt,phi,eps1,t1,AA,nn,G01)
    lamda1=2.*mu1.*v./(1-2.*v);  
    C1=trapz(r,(2.*pi.*r.*(phi(1,:).^2).*(lamda1+2.*mu1)));
    k1=Ep.*Ap+trapz(r,2.*pi.*r.*(lamda1+2.*mu1).*(phi(1,:).^2));
    m1=trapz(r,(-2.*pi.*r.*mu1.*(phi(2,:).^2)));
    s1=-Ep.*Ap-trapz(r,(2.*pi.*r.*(lamda1+2.*mu1).*(phi(1,:).^2)));
    f1=0; 
   mu2=Gs12(z11,z2,rt,phi,eps2,t2,AA,nn,G02);
    lamda2=2.*mu2.*v./(1-2.*v);
       
    C2=trapz(r,(2.*pi.*r.*(phi(1,:).^2).*(lamda2+2.*mu2)));
    k2=Ep.*Ap+trapz(r,2.*pi.*r.*(lamda2+2.*mu2).*(phi(1,:).^2));
    m2=trapz(r,(-2.*pi.*r.*mu2.*(phi(2,:).^2)));
    s2=-Ep.*Ap-trapz(r,(2.*pi.*r.*(lamda2+2.*mu2).*(phi(1,:).^2)));
    f2=0;     

C3=0;k3=0;m3=0;s3=0;f3=0;z0=0;

    mu3=Gs13(z22,z3,rt,phi,eps3,t3,AA,nn,G03);
    lamda3=2.*mu3.*v./(1-2.*v);
       
    C3=trapz(r,(2.*pi.*r.*(phi(1,:).^2).*(lamda3+2.*mu3)));
    k3=Ep.*Ap+trapz(r,2.*pi.*r.*(lamda3+2.*mu3).*(phi(1,:).^2));
    m3=trapz(r,(-2.*pi.*r.*mu3.*(phi(2,:).^2)));
    s3=-Ep.*Ap-trapz(r,(2.*pi.*r.*(lamda3+2.*mu3).*(phi(1,:).^2)));
    f3=0; 

   C4=0;k4=0;m4=0;s4=0;f4=0;z0=0;

    mu4=Gs14(z33,z4,rt,phi,eps4,t4,AA,nn,G04);
    lamda4=2.*mu4.*v./(1-2.*v);  
    C4=trapz(r,(2.*pi.*r.*(phi(1,:).^2).*(lamda4+2.*mu4)));
    k4=Ep.*Ap+trapz(r,2.*pi.*r.*(lamda4+2.*mu4).*(phi(1,:).^2));
    m4=trapz(r,(-2.*pi.*r.*mu4.*(phi(2,:).^2)));
    s4=-Ep.*Ap-trapz(r,(2.*pi.*r.*(lamda4+2.*mu4).*(phi(1,:).^2)));
    f4=0; 
     
%phi coefficients 

ms11=zeros(1,N);ms12=zeros(1,N);

for i=1:N

       mu1=Gs21(rt(i),rt,phi,eps1,t1,AA,nn,G01);
       lamda1=2.*mu1.*v./(1-2.*v);
       ms11(i)=trapz(z1,(1+(r.*mu1./mu1)));
       ms12(i)=-trapz(z1,((lamda1+2.*mu1).*(t1.^2)))./trapz(z1,(mu1.*...
           (eps1(1,:).^2)));
end

ms21=zeros(1,N);ms22=zeros(1,N);
for i=1:N
           mu2=Gs22(rt(i),rt,phi,eps2,t2,AA,nn,G02);
       lamda2=2.*mu2.*v./(1-2.*v);
       ms21(i)=trapz(z2,(1+(r.*mu2./mu2)));
       ms22(i)=-trapz(z2,((lamda2+2.*mu2).*(t2.^2)))./trapz(z2,(mu2.*...
           (eps2(1,:).^2)));
end

ms31=zeros(1,N);ms32=zeros(1,N);
for i=1:N
       mu3=Gs23(rt(i),rt,phi,eps3,t3,AA,nn,G03);
       lamda3=2.*mu3.*v./(1-2.*v);
      
       ms31(i)=trapz(z3,(1+(r.*mu3./mu3)));
       ms32(i)=-trapz(z3,((lamda3+2.*mu3).*(t3.^2)))./trapz(z3,(mu3.*...
           (eps3(1,:).^2)));
end

ms41=zeros(1,N);ms42=zeros(1,N);
for i=1:N
       mu4=Gs24(rt(i),rt,phi,eps4,t4,AA,nn,G04);
       lamda4=2.*mu4.*v./(1-2.*v);
       ms41(i)=trapz(z4,(1+(r.*mu4./mu4)));
       ms42(i)=-trapz(z4,((lamda4+2.*mu4).*(t4.^2)))./trapz(z4,(mu4.*...
           (eps4(1,:).^2)));
end
ms1=zeros(1,N);ms2=zeros(1,N);ms3=zeros(1,N);ms4=zeros(1,N);
ms1=ms11+ms21+ms31+ms41; 
ms2=ms12+ms22+ms32+ms42;
B=zeros(n4,n4);
B=zeros(n1,n1);
B(1,1)=(2.*f1.*k1./(dz.*s1))-(2.*k1./(dz.^2))-(C1.*f1./s1)+m1; %equ1
B(1,2)=2.*k1./(dz.^2);
B(n1,n1)=1; %equ6
B(n1,n1+1)=-1;
%--------------------------------------------------------------------------
B(n2,n2)=1; %equ11
B(n2,n2+1)=-1;

%--------------------------------------------------------------------------
B(n3,n3)=1;%equ16
B(n3,n3+1)=-1;
% %------------------------------------------------------------------------
B(n4,n4)=1;%equ21
%--------------------------------------------------------------------------

for i=2:n1-1
    for j=i-1:i-1
    B(i,j)=k1./(dz.^2)-C1./(2.*dz);
    end
end
 for i=2:n1-1
    for j=i:i
    B(i,j)=-2*k1./(dz.^2)+m1;
end
end

for i=2:n1-1
   for j=i+1:i+1
    B(i,j)=k1./(dz.^2)+C1./(2.*dz);
   end
end

%--------------------------------------------------------------------------
 for i=n1+1:n2-1
    for j=i-1:i-1
  B(i,j)=k2./(dz.^2)-C2./(2.*dz);
end
end

 for i=n1+1:n2-1
  for j=i:i
    B(i,j)=-2*k2./(dz.^2)+m2;
end
end

for i=n1+1:n2-1
   for j=i+1:i+1
      B(i,j)=k2./(dz.^2)+C2./(2.*dz);
end
end
%--------------------------------------------------------------------------
for i=n2+1:n3-1
    for j=i-1:i-1
      B(i,j)=k3./(dz.^2)-C3./(2.*dz);
    end
end
for i=n2+1:n3-1
  for j=i:i
     B(i,j)=-2*k3./(dz.^2)+m3;
end
end

for i=n2+1:n3-1
    for j=i+1:i+1
    B(i,j)=k3./(dz.^2)+C3./(2.*dz);
end
end

%--------------------------------------------------------------------------
for i=n3+1:n4-1
    for j=i-1:i-1
      B(i,j)=k4./(dz.^2)-C4./(2.*dz);
    end
end
for i=n3+1:n4-1
    for j=i:i
     B(i,j)=-2*k4./(dz.^2)+m4;
end
end

for i=n3+1:n4-1
    for j=i+1:i+1
  B(i,j)=k4./(dz.^2)+C4./(2.*dz);
    end 
end

%--------------------------------------------------------------------------
% B=B(1:n2,1:n2)
B=B(1:n4,1:n4);
f(1,1)=(P0.*C1./s1)-2.*P0.*k1./(s1.*dz);
f(2:n4)=0;
u=inv(B)*f'
eps1(1,:)=u(1:n1,1);
eps2(1,:)=u(n1+1:n2,1);
eps3(1,:)=u(n2+1:n3,1);
eps4(1,:)=u(n3+1:n4,1); 
t1=dfdx(eps1(1,:),z1,N);
t2=dfdx(eps2(1,:),z2,N);
t3=dfdx(eps3(1,:),z3,N);
t4=dfdx(eps4(1,:),z4,N);
y1new=ms1;
y2new=ms2;
y1old=y1;y2old=y2;
my=(y1old-y1new)./y1old;
end
dlmwrite('axial1.txt',mu1)
dlmwrite('axial2.txt',mu2);
dlmwrite('axial3.txt',mu3);
dlmwrite('axial4.txt',mu4);
end

%--------------------------------------------------------------------------
function dphidr =phiode(r,phi,rt,y1,y2)
x1=interp1(rt,y1,r,'cubic');
x2=interp1(rt,y2,r,'cubic');
  dphidr = [ 
     
           phi(2) 
           -phi(2).*(x1.*1./r)-phi(1).*x2
            ];
end

function res = phibc(phia,phib)
  res = [ phia(1)-1 
           phib(1)

          ];
end

 %-------------------------------------------------------------------------
function finddmdr=dfdx(mm,rt,N)
        diff1=(mm(1)-mm(2))./(rt(1)-rt(2))+(mm(1)-mm(3))./(rt(1)-rt(3))-...
            (mm(2)-mm(3))./(rt(2)-rt(3));
        diffi=(mm(1:(N-2))-mm(2:(N-1)))./(rt(1:(N-2))-rt(2:(N-1)))....
            -(mm(1:(N-2))-mm(3:N))./(rt(1:(N-2))-rt(3:N))....
            +(mm(2:(N-1))-mm(3:N))./(rt(2:(N-1))-rt(3:N));
        diffn=-(mm(N-1)-mm(N-2))./(rt(N-1)-rt(N-2))+(mm(N-2)-...
            mm(N))./(rt(N-2)-rt(N))+(mm(N-1)-mm(N))./(rt(N-1)-rt(N));
finddmdr=[diff1,diffi,diffn];
end
%--------------------------------------------------------------------------
function findG11=Gs11(z0,z1,rt,phi,eps1,t1,AA,nn,G01);
epsiz1=interp1(z1,eps1(1,:),z0);
depsiz1=interp1(z1,t1,z0);
epsilon_r1=0;epsilon_theta1=0;epsilon_rtheta1=0;epsilon_thetaz1=0;
epsilon_z1=depsiz1.*phi(1,:);
epsilon_zr1=phi(2,:).*epsiz1;
epsilon_q1=2/3.*sqrt(((epsilon_z1).^2+epsilon_theta1.^2+...
    epsilon_r1.^2+epsilon_z1.^2)+3.*(epsilon_rtheta1.^2+epsilon_thetaz1.^2+...
    epsilon_zr1.^2));
findG11=AA.*epsilon_q1.^nn;
mm=find(findG11>G01);
findG11(mm)=G01;
end

function findG12=Gs12(z11,z2,rt,phi,t2,eps2,AA,nn,G02)
epsiz2=interp1(z2,eps2(1,:),z11);
depsiz2=interp1(z2,t2,z11);
epsilon_r2=0;epsilon_theta2=0;epsilon_rtheta2=0;epsilon_thetaz2=0;

epsilon_z2=depsiz2.*phi(1,:);


epsilon_zr2=phi(2,:).*epsiz2;


epsilon_q2=2/3.*sqrt(((epsilon_z2).^2+epsilon_theta2.^2+...
    epsilon_r2.^2+epsilon_z2.^2)+3.*(epsilon_rtheta2.^2+epsilon_thetaz2.^2+...
    epsilon_zr2.^2));

findG12=AA.*epsilon_q2.^nn;
mm=find(findG12>G02);
findG12(mm)=G02;

end

function findG13=Gs13(z22,z3,rt,phi,t3,eps3,AA,nn,G03)

epsiz3=interp1(z3,eps3(1,:),z22);
depsiz3=interp1(z3,t3,z22);
epsilon_r3=0;epsilon_theta3=0;epsilon_rtheta3=0;epsilon_thetaz3=0;
epsilon_z3=depsiz3.*phi(1,:);
epsilon_zr3=phi(2,:).*epsiz3;
epsilon_q3=2/3.*sqrt(((epsilon_z3).^2+epsilon_theta3.^2+...
    epsilon_r3.^2+epsilon_z3.^2)+3.*(epsilon_rtheta3.^2+epsilon_thetaz3.^2+...
    epsilon_zr3.^2));
findG13=AA.*epsilon_q3.^nn;
mm=find(findG13>G03);
findG13(mm)=G03;

end

function findG14=Gs14(z33,z4,rt,phi,t4,eps4,AA,nn,G04)
epsiz4=interp1(z4,eps4(1,:),z33);
depsiz4=interp1(z4,t4,z33);
epsilon_r4=0;epsilon_theta4=0;epsilon_rtheta4=0;epsilon_thetaz4=0;
epsilon_z4=depsiz4.*phi(1,:);
epsilon_zr4=phi(2,:).*epsiz4;
epsilon_q4=2/3.*sqrt(((epsilon_z4).^2+epsilon_theta4.^2+...
    epsilon_r4.^2+epsilon_z4.^2)+3.*(epsilon_rtheta4.^2+epsilon_thetaz4.^2+...
    epsilon_zr4.^2));

findG14=AA.*epsilon_q4.^nn;
mm=find(findG14>G04);
findG14(mm)=G04;
end
%--------------------------------------------------------------------------

function findG21=Gs21(rp,rt,phi,t1,eps1,AA,nn,G01)
phi_r=interp1(rt,phi(1,:),rp);
dphi_r=interp1(rt,phi(2,:),rp);
epsilon_r1=0;epsilon_theta1=0;epsilon_rtheta1=0;epsilon_thetaz1=0;
epsilon_z1=t1.*phi_r;
epsilon_zr1=dphi_r.*eps1(1,:);
epsilon_q1=2/3.*sqrt(1/2.*((epsilon_r1-epsilon_theta1).^2+(epsilon_z1-...
epsilon_theta1).^2+(epsilon_r1-epsilon_z1).^2+epsilon_theta1.^2+...
epsilon_r1.^2+epsilon_z1.^2)+3.*(epsilon_rtheta1.^2+epsilon_thetaz1.^2+...
epsilon_zr1.^2));
findG21=AA.*epsilon_q1.^nn;
mm=find(findG21>G01);
findG21(mm)=G01;
end

function findG22=Gs22(rp,rt,phi,t2,eps2,AA,nn,G02)

phi_r=interp1(rt,phi(1,:),rp);
dphi_r=interp1(rt,phi(2,:),rp);
epsilon_r2=0;epsilon_theta2=0;epsilon_rtheta2=0;epsilon_thetaz2=0;
epsilon_z2=t2.*phi_r;
epsilon_zr2=dphi_r.*eps2(1,:);
epsilon_q2=2/3.*sqrt(1/2.*((epsilon_r2-epsilon_theta2).^2+(epsilon_z2-...
epsilon_theta2).^2+(epsilon_r2-epsilon_z2).^2+epsilon_theta2.^2+...
epsilon_r2.^2+epsilon_z2.^2)+3.*(epsilon_rtheta2.^2+epsilon_thetaz2.^2+...
epsilon_zr2.^2));
findG22=AA.*epsilon_q2.^nn;
mm=find(findG22>G02);
findG22(mm)=G02;
end

function findG23=Gs23(rp,rt,phi,t3,eps3,AA,nn,G03)

phi_r=interp1(rt,phi(1,:),rp);
dphi_r=interp1(rt,phi(2,:),rp);
epsilon_r3=0;epsilon_theta3=0;epsilon_rtheta3=0;epsilon_thetaz3=0;
epsilon_z3=t3.*phi_r;
epsilon_zr3=dphi_r.*eps3(1,:);
epsilon_q3=2/3.*sqrt(1/2.*((epsilon_r3-epsilon_theta3).^2+(epsilon_z3-...
epsilon_theta3).^2+(epsilon_r3-epsilon_z3).^2+epsilon_theta3.^2+...
epsilon_r3.^2+epsilon_z3.^2)+3.*(epsilon_rtheta3.^2+epsilon_thetaz3.^2+...
epsilon_zr3.^2));
findG23=AA.*epsilon_q3.^nn;
mm=find(findG23>G03);
findG23(mm)=G03;

end

function findG24=Gs24(rp,rt,phi,t4,eps4,AA,nn,G04)
phi_r=interp1(rt,phi(1,:),rp);
dphi_r=interp1(rt,phi(2,:),rp);
epsilon_r4=0;epsilon_theta4=0;epsilon_rtheta4=0;epsilon_thetaz4=0;
epsilon_z4=t4.*phi_r;
epsilon_zr4=dphi_r.*eps4(1,:);
epsilon_q4=2/3.*sqrt(1/2.*((epsilon_r4-epsilon_theta4).^2+(epsilon_z4-...
epsilon_theta4).^2+(epsilon_r4-epsilon_z4).^2+epsilon_theta4.^2+...
epsilon_r4.^2+epsilon_z4.^2)+3.*(epsilon_rtheta4.^2+epsilon_thetaz4.^2+...
epsilon_zr4.^2));
findG24=AA.*epsilon_q4.^nn;
mm=find(findG24>G04);
findG24(mm)=G04;


end