How to surf plot and line plot together for the code?

Question

0 votes

I need to plot like

But for my code it showing some error

[x, y] = meshgrid(-.5:0.01:.5);
m =1;
sgm =1;
wo = 0.01;
lambda=532 *10^-9;
k=2*pi/lambda;
 
Omega = .2;
%OTOPS model
z=20 
Pho=0.007 
 for u=1:length(z)
    I=0;
    Gamma = ((1i.*k)./(2.*z(u))) + (1./(Pho(u).^2)) + (1./(wo.^2));
    Gamma_star = subs(Gamma, 1i, -1i);
    alpha = (Gamma_star) - 1./((Gamma.*Pho(u).^4));
    C = (1./(4.*(lambda.^2).*(z(u).^2))).*((pi.^2)./(alpha.*Gamma)).*(1./(2.*1i.*sqrt(Gamma)).^m); 
    fx = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
    fy = ((Omega.^2)./(4.*Gamma)) + ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
    gx = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*x.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*x.^2)./(4.*z(u).^2.*Gamma));
    gy = ((Omega.^2)./(4.*Gamma)) - ((1i.*k.*y.*Omega)./(2.*z(u).*Gamma)) - ((k.^2.*y.^2)./(4.*z(u).^2.*Gamma));
    Bx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    Bx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    By_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    By_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) + (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    Kx_p = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    Kx_n = ((-1i.*k.*x)./(2.*z(u))) + ((1i.*k.*x)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    Ky_p = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) + (Omega./2);
    Ky_n = ((-1i.*k.*y)./(2.*z(u))) + ((1i.*k.*y)./(2.*z(u).*Gamma.*Pho(u).^2)) - (Omega)./(2.*Gamma.*Pho(u).^2) - (Omega./2);
    
    X = 0;
    for l = 0:m
      L = (((1i).*sgm).^(m-l)).*nchoosek(m,l);
      sum_j = 0;
      for j = 0:m
          J= (((-1i).*sgm).^(m-j)).*nchoosek(m,j);  
          
          sum1 = 0.0;
          for p = 0:1/2:l/2
                faktor1_P = (((-1).^p).*factorial(l))./(gamma(p+1).*factorial(l-(2.*p))).*(((2.*1i)./(sqrt(Gamma))).^(l-(2.*p))).*exp(fx).*exp(fy);
                for s1 = 0:l-2*p
                    faktor1_s1 =  nchoosek((l-(2.*p)),s1).*((1./(Pho(u).^2)).^s1).*((((1i.*k.*x)./(2.*z(u)))+(Omega./2)).^(l-(2.*p)-s1)).*(1./(((2.*1i.*sqrt(alpha))).^(j+s1)));
                    for q = 0:1/2:(m-l)/2
                        faktor1_Q = (((-1).^q).*factorial(m-l))./(gamma(q+1).*factorial(m-l-(2.*q))).*(((2.*1i)./(sqrt(Gamma))).^(m-l-(2.*q)));
                        sum_inner_1 = 0;
                        for s2 = 0:m-l-2*q
                            E1 = (exp(((Bx_p).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_p)./(sqrt(alpha)))).*exp(((By_p).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_p)./(sqrt(alpha))))) - (exp(((Bx_n).^2)./alpha).*hermiteH(j+s1,((1i.*Bx_n)./(sqrt(alpha)))).*exp(((By_n).^2)./alpha).*hermiteH(m-j+s2,((1i.*By_n)./(sqrt(alpha)))));
                            sum_inner_1 = sum_inner_1 + nchoosek((m-l-(2.*q)),s2).*((1./(Pho(u).^2)).^s2).*((((1i.*k.*y)./(2.*z(u)))+(Omega./2)).^(m-l-(2.*q)-s2)).*(1./(((2.*1i.*sqrt(alpha))).^(m-j+s2))).*E1;
                        end
                    end
                end
                sum1 = sum1 + faktor1_P.*faktor1_s1.*faktor1_Q.*sum_inner_1;
           end
      
 
       sum_j= sum_j + J.*(sum1);
      end
     X = X + (L.*sum_j);
    end  
    
    X  = (C.*X);  
    I = I + abs(vpa(X)) ;
    figure;
    surf(x, y, I);
    hold on; 
    plot3(x(1,:), y(1,:), I(51,:), 'r-', 'LineWidth', 2);
    plot3(y(1,:), x(1,:), I(51,:), 'r-', 'LineWidth', 2);
    hold off;  
    view(-45, 30);
    xlim([-.5, .5]);
    ylim([-.5, .5]);
    legend('Surface', 'Line');
end