clc
clear all;
d=0.046;C_1=0.002;t=0.4;p_1=7870;p_2=865;u=0.00003;b=0.0001;l=5;A_2=1;
w=0:0.5:25;
z_n=0:2.5/75:2.5;
Q=zeros(51,76);
for i=1:51
for j=1:76
Q1(i,j)=(2*b.*exp((3.*z_n(j))./(2.*l) - (t.*w(i))./2))/(w(i)+ (6.*d.^3.*u)./(C_1.^3.*((3.*p_2)./5000 + p_1.*(d - 3./100).^2))) - (2.*b.*w(i).*exp((3.*z_n(j))./(2.*l) - (t.*w(i))./2))./(w(i)+ (6.*d.^3.*u)./(C_1.^3.*((3.*p_2)./5000 + p_1.*(d - 3./100).^2))).^2 - (b.*t.*w(i).*exp((3.*z_n(j))./(2.*l) - (t.*w(i))./2))/(w(i) + (6.*d.^3.*u)./(C_1.^3.*((3.*p_2)./5000 + p_1.*(d - 3./100).^2)));
Q2(i,j)=(3.*A_2.*b.*w(i).*exp((3.*z_n(j))./(2.*l) - (t.*w(i))./2))./(l.*(w(i)+ (6.*d.^3.*u)./(C_1.^3.*((3.*p_2)./5000 + p_1.*(d - 3./100).^2))));
end
end
[z_nz_nw,www]=meshgrid(z_n,w);
subplot(1,2,1)
surf(Q1);
view(50,40)
grid on
xlabel('z_n');
ylabel('w');
zlabel('a/A_2');
hold on;
surf(Q2);
view(50,40)
grid on
xlabel('z_n');
ylabel('w');
zlabel('a/A_2');
samedot=(abs(Q1(i,j)-Q2(i,j))<=0.5);
z_nz_n=samedot.*z_n;
ww=samedot.*w;
zz=samedot.*Q1(i.j);
plot3(z_nz_n(samedot~=0),ww(samedot~=0),zz(samedot~=0),'r-');
subplot(1,2,2)
h1=plot3(z_nz_n(samedot~=0),ww(samedot~=0),zz(samedot~=0),'r-');