my code is not running , i have attached my code.

%In image, there are two equations. I put equation (3) {Zc} in equation (5) and get quadratic equation of ZL . I developed a code for solution of quadratic equation , in which i used 'if ' condition {for eq. 3 ,where E>0 rho _w =0) .Please look at it and correct it. Result is not getting ,what it is expected,may be something is wrong.

%thanks in Advance

%link for data

% https://drive.google.com/drive/folders/1fmQHWyp_cvBs9ULuWUpL3uCuSAykJbZl?usp=sharing

E=xlsread('toplatlongecs','c1:c292681');
N=xlsread('geoidlatlongecs','c1:c292681');

% syms N E real
L0=2400;  %load in metres
r_c=2670; % density of crust in kg/m3
r_w=1030; % density of water in kg/m3
r_a=3200; % density of asthenosphere in kg/m3
r_m=3300; % density of lithosphere mantle in kg/m3
Zmax=300000; %compensation level, zmax in metres
g=9.8;       % gravity in m/s2
r_ac=r_a-r_c; % density contrast b/w asthenosphere and crust
r_cw=r_c-r_w; % density contrast b/w  crust and water
r_wc=r_w-r_c; % density contrast b/w water and crust
r_ma=r_m-r_a; % density contrast b/w mantle lithosphere and asthenosphere
r_mc=r_m-r_c; % density contrast b/w mantle lithosphere  and crust
rmc=(r_mc.^2);
r_cm=r_c-r_m; 
rma=(r_ma.^2);
ur_ac=1/(r_ac);
ur_mc=1/(r_mc);
rcw=(r_cw.^2);
ZL= 120000; %in metres
ZC= 27000;   % in metres
Zmax2=Zmax^2;
E2=E.^2;
ZC2=ZC.^2;
ZL2=ZL.^2;
C2=3.14*6.67*10^(-11);
N0=-(C2*(E2*r_w+ (ZC2-E2)*r_c+(ZL2-ZC2)*r_m+ (Zmax2-ZL2)*r_a))./ 9.8 -N ;
N22=-(N0+N)*9.8/C2;

%r_cw=r_c;
%zC=(r_a*L0+E*(r_cw+zL*r_ma)*ur_mc;

%A=(r_c*rma/rmc)+r_m-(r_m*rma/rmc)-r_a;
for i= 1:292681;
if E(i)>0
    
A=r_ma+(rma*r_cm)/(rmc);
B(i)=(2*r_ma/rmc)*(E(i)*r_c+L0*r_a)*(r_cm);
D1(i)=(E(i).^2)*(r_wc+(r_cm*r_c*r_c/rmc));
D2(i)=(2*r_a*r_c*r_cm*L0*E(i))/(rmc);
D3=Zmax2*r_a;
D4=(r_a*r_a*L0*L0*r_cm)/rmc;   
D5(i)=D1(i)+D2(i)+D3+D4;
A2=2*A;

 N23(i)=N22(i)';  
 c11(i)=D5(i)-N23(i);
 xL11(i)=(-(B(i))+((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
 %xL1=real(xL11);
 xL1(i)=xL11(i)/1000;
 xL33(i)=xL1(i);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 xL22(i)=(-(B(i))-((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
% xL2=real(xL22);
 xL2(i)=xL22(i)/1000;
 xL44(i)=xL2(i);
       
else
A=r_ma+(rma*r_cm)/(rmc);
B(i)=(2*r_ma/rmc)*(E(i)*r_cw+L0*r_a)*(r_cm);
D1(i)=(E(i).^2)*(r_wc+(r_cm*rcw/rmc));
D2(i)=(2*r_a*r_cw*r_cm*L0*E(i))/(rmc);
D3=Zmax2*r_a;
D4=(r_a*r_a*L0*L0*r_cm)/rmc;
D5(i)=D1(i)+D2(i)+D3+D4;
A2=2*A;



 N23(i)=N22(i)';  
 c11(i)=D5(i)-N23(i);
 xL11(i)=(-(B(i))+((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
 %xL1=real(xL11);
 xL1(i)=xL11(i)/1000;
 xL33(i)=xL1(i);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 xL22(i)=(-(B(i))-((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
% xL2=real(xL22);
 xL2(i)=xL22(i)/1000;
 xL44(i)=xL2(i);
 
end 
end;
 
% xL11=(-(B)+((B.^2)-4*A*c11).^(1/2))/(A2);
%  %xL1=real(xL11);
%  xL1=xL11;
%  
%  xL33=xL1/1000;
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  xL22=(-(B)-((B.^2)-4*A*c11).^(1/2))/(A2);
% % xL2=real(xL22);
%  xL2=xL22/1000;
%  
%  xL44=xL2;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 zL2=reshape(xL33,541,541);
 zL=zL2';
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 xL334=xL33';
 zC2=(r_a*L0+E*r_c+xL334*r_ma)/(r_mc);
 zC3=reshape(zC2,541,541);
 zC=zC3';
 
 subplot(2,1,1)
im= imagesc([117 135], [-38 -20] ,zL);
 colorbar
 hold on          

 colorbar; 

%brighten('h',);
%C=contour(A,'LineColor','black');
 x=117:0.033333333:135;
 y=-38:0.033333333:-20;
 
 minx = round(min(zL(:)),1);
 maxx = round(max(zL(:)),1);
 levels =  minx:5:maxx;
 levels2=round(levels);
 
[X,Y] = meshgrid(x,y);
%[C,h]=contour(X,Y,csia21,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
[C,h]=contour(X,Y,zL,levels2,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
%C=contour(A,'ShowText','on','LineColor','black');
clabel(C,h,'FontSize',6);

subplot(2,1,2)
%jet map 
im= imagesc([117 135], [-38 -20] ,zC);
 colorbar
 hold on          
colorbar; 

%brighten('h',);
%C=contour(A,'LineColor','black');
 x=117:0.033333333:135;
 y=-38:0.033333333:-20;
 
 minx = round(min(zC(:)),1);
 maxx = round(max(zC(:)),1);
 levels =  minx:5:maxx;
 levels2=round(levels);
 
[X,Y] = meshgrid(x,y);
%[C,h]=contour(X,Y,csia21,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
[C,h]=contour(X,Y,zC,levels2,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
%C=contour(A,'ShowText','on','LineColor','black');
clabel(C,h,'FontSize',6);

[Image Analyst]

1. Can you attach your data here with the paperclip icon? Or is it too big (more than 5 MB)?
2. Can you format your code as code?
3. Can you describe what is wrong. I imagine that I will run it and it will run and produce some numbers but that you don't like them for some unspecified reason.
4. What is the expected answer?

Please read this link:

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[Priyank Pathak]

@Image Analyst,data is more than 5MB. That's why i attached a link.

% i attached an image of required result. left one for zC and right one for zL.

%link for data: LINK

% not getting required result.

E=xlsread('toplatlongecs','c1:c292681');
N=xlsread('geoidlatlongecs','c1:c292681');

% syms N E real
L0=2400;  %load in metres
r_c=2670; % density of crust in kg/m3
r_w=1030; % density of water in kg/m3
r_a=3200; % density of asthenosphere in kg/m3
r_m=3300; % density of lithosphere mantle in kg/m3
Zmax=300000; %compensation level, zmax in metres
g=9.8;       % gravity in m/s2
r_ac=r_a-r_c; % density contrast b/w asthenosphere and crust
r_cw=r_c-r_w; % density contrast b/w  crust and water
r_wc=r_w-r_c; % density contrast b/w water and crust
r_ma=r_m-r_a; % density contrast b/w mantle lithosphere and asthenosphere
r_mc=r_m-r_c; % density contrast b/w mantle lithosphere  and crust
rmc=(r_mc.^2);
r_cm=r_c-r_m; 
rma=(r_ma.^2);
ur_ac=1/(r_ac);
ur_mc=1/(r_mc);
rcw=(r_cw.^2);
ZL= 120000; %in metres
ZC= 27000;   % in metres
Zmax2=Zmax^2;
E2=E.^2;
ZC2=ZC.^2;
ZL2=ZL.^2;
C2=3.14*6.67*10^(-11);
N0=-(C2*(E2*r_w+ (ZC2-E2)*r_c+(ZL2-ZC2)*r_m+ (Zmax2-ZL2)*r_a))./ 9.8 -N ;
N22=-(N0+N)*9.8/C2;


for i= 1:292681;
if E(i)>0
    
A=r_ma+(rma*r_cm)/(rmc);
B(i)=(2*r_ma/rmc)*(E(i)*r_c+L0*r_a)*(r_cm);
D1(i)=(E(i).^2)*(r_wc+(r_cm*r_c*r_c/rmc));
D2(i)=(2*r_a*r_c*r_cm*L0*E(i))/(rmc);
D3=Zmax2*r_a;
D4=(r_a*r_a*L0*L0*r_cm)/rmc;   
D5(i)=D1(i)+D2(i)+D3+D4;
A2=2*A;

 N23(i)=N22(i)';  
 c11(i)=D5(i)-N23(i);
 xL11(i)=(-(B(i))+((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
 %xL1=real(xL11);
 xL1(i)=xL11(i)/1000;
 xL33(i)=xL1(i);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 xL22(i)=(-(B(i))-((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
 % xL2=real(xL22);
 xL2(i)=xL22(i)/1000;
 xL44(i)=xL2(i);
       
else
A=r_ma+(rma*r_cm)/(rmc);
B(i)=(2*r_ma/rmc)*(E(i)*r_cw+L0*r_a)*(r_cm);
D1(i)=(E(i).^2)*(r_wc+(r_cm*rcw/rmc));
D2(i)=(2*r_a*r_cw*r_cm*L0*E(i))/(rmc);
D3=Zmax2*r_a;
D4=(r_a*r_a*L0*L0*r_cm)/rmc;
D5(i)=D1(i)+D2(i)+D3+D4;
A2=2*A;

N23(i)=N22(i)';  
 c11(i)=D5(i)-N23(i);
 xL11(i)=(-(B(i))+((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
 %xL1=real(xL11);
 xL1(i)=xL11(i)/1000;
 xL33(i)=xL1(i);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 xL22(i)=(-(B(i))-((B(i).^2)-4*A*c11(i)).^(1/2))/(A2);
% xL2=real(xL22);
 xL2(i)=xL22(i)/1000;
 xL44(i)=xL2(i);
 
end 
end;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 zL2=reshape(xL33,541,541);
 zL=zL2';
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 xL334=xL33';
 zC2=(r_a*L0+E*r_c+xL334*r_ma)/(r_mc);
 zC3=reshape(zC2,541,541);
 zC=zC3';
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,1)
im= imagesc([117 135], [-38 -20] ,zL);
 colorbar
 hold on          

 colorbar; 

 x=117:0.033333333:135;
 y=-38:0.033333333:-20;
 
 minx = round(min(zL(:)),1);
 maxx = round(max(zL(:)),1);
 levels =  minx:5:maxx;
 levels2=round(levels);
 
[X,Y] = meshgrid(x,y);
[C,h]=contour(X,Y,zL,levels2,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
%C=contour(A,'ShowText','on','LineColor','black');
clabel(C,h,'FontSize',6);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,1,2) 
im= imagesc([117 135], [-38 -20] ,zC);
 colorbar
 hold on          
colorbar; 

 x=117:0.033333333:135;
 y=-38:0.033333333:-20;
 
 minx = round(min(zC(:)),1);
 maxx = round(max(zC(:)),1);
 levels =  minx:5:maxx;
 levels2=round(levels);
 
[X,Y] = meshgrid(x,y);
[C,h]=contour(X,Y,zC,levels2,'linewidth',0.2,'LineColor',[0,0,0]+0.5);
clabel(C,h,'FontSize',6);