How to resolve the Plotting issue in my code?

I done the partial differentiation of W and W3 and then substitute with an array of values and tried to plot.

clc;close all; clear all;
w=-3;T=38;S=28;Om =1;  w0 =0.002; sigma0 = 0.04;muu0=4*pi*10^-7;
lambda=532*10^-9;k=2*pi/lambda;z=100;epsilon0=8.85*10^-12;
alphac=0.02;T=3.609699516774368e+03;
syms  r1x r1y r2x r2y 
J = 1i*k/(2*z) - T;
J1 = conj(J) ;
A = 1/(2*w0) + 1/(2*sigma0^2) + T;
Bx = (1i*k/2*z)*(r1x+r2x) - T*(r1x - r2x);By = (1i*k/2*z)*(r1y+r2y) - T*(r1y - r2y);
C = 2/(w0^2) + k^2/(4*A*z^2);
Dx = (1i*k/z)*(r1x-r2x) + 2*Om; Dy = (1i*k/z)*(r1y-r2y) + 2*Om;
Dx1 = (1i*k/z)*(r1x-r2x) - 2*Om; Dy1 = (1i*k/z)*(r1y-r2y) - 2*Om;
Ex = 0.5*(Dx - (1i*k*Bx)/(2*A*z)); Ey = 0.5*(Dy - (1i*k*By)/(2*A*z));
Ex1 = 0.5*(Dx1 - (1i*k*Bx)/(2*A*z)); Ey1 = 0.5*(Dy1 - (1i*k*By)/(2*A*z));
Fx=Bx + Om;Fy = By + Om;
Fx1=Bx - Om;Fy1 = By - Om;
Gx=0.5*((1i*k/z)*(r1x-r2x) - (1i*k*Fx/2*A*z) );Gy=0.5*((1i*k/z)*(r1y-r2y) - (1i*k*Fy/2*A*z) );
Gx1=0.5*((1i*k/z)*(r1x-r2x) - (1i*k*Fx1/2*A*z) );Gy1=0.5*((1i*k/z)*(r1y-r2y) - (1i*k*Fy1/2*A*z) );
H = A + (k^2*w0^2)/(8*z^2);
Ix = 0.5*(Bx - (1i*k*w0*w0*Dx)/(4*z));Iy = 0.5*(By - (1i*k*w0*w0*Dy)/(4*z));
Ix1 = 0.5*(Bx - (1i*k*w0*w0*Dx1)/(4*z));Iy1 = 0.5*(By - (1i*k*w0*w0*Dy1)/(4*z));
Jx = 0.5*(Fx + (k*k*w0*w0*(r1x-r2x))/(4*z*z));Jy = 0.5*(Fy + (k*k*w0*w0*(r1y-r2y))/(4*z*z));
Jx1 = 0.5*(Fx1 + (k*k*w0*w0*(r1x-r2x))/(4*z*z));Jy1 = 0.5*(Fy1 + (k*k*w0*w0*(r1y-r2y))/(4*z*z));
M1 =(( Ex^2 + Ey^2 )/(C^2) + (1/C) - (Ix^2 + Iy^2 )/(4*H^2) - (1/4*H) + (1i*(Ix*Ey - Ex*Iy)/(C*H)))*exp(((Bx^2 + By^2)/(4*A)) + ((Ex^2 + Ey^2)/C));
M2 =(( Ex1^2 + Ey1^2 )/(C^2) + (1/C) - (Ix1^2 + Iy1^2 )/(4*H^2) - (1/4*H) + (1i*(Ix1*Ey1 - Ex1*Iy1)/(C*H)))*exp(((Bx^2 + By^2)/(4*A)) + ((Ex1^2 + Ey1^2)/C));
M3 =(( Gx^2 + Gy^2 )/(C^2) + (1/C) - (Jx^2 + Jy^2 )/(4*H^2) - (1/4*H) + (1i*(Jx*Gy - Gx*Jy)/(C*H)))*exp(((Fx^2 + Fy^2)/(4*A)) + ((Gx^2 + Gy^2)/C));
M4 =(( Gx1^2 + Gy1^2 )/(C^2) + (1/C) - (Jx1^2 + Jy1^2 )/(4*H^2) - (1/4*H) + (1i*(Jx1*Gy1 - Gx1*Jy1)/(C*H)))*exp(((Fx1^2 + Fy1^2)/(4*A)) + ((Gx1^2 + Gy1^2)/C));
P = (k*k)/(16*A*C*z*z);
W=P*exp(conj(J)*(r1x^2+r1y^2) + (J)*(r2x^2+r2y^2) + 2*T*(r1x*r2x + r1y*r2y))*(M1 + M2 - M3 - M4);
dW1=diff(W,r2x);
dW2=diff(W,r2y);
W3 = r1y*dW1 - r1x*dW2;
r=linspace(-1,1,100);
q1 =(subs(W3,{r1x,r1y,r2x,r2y},{r,r,r,r}));
q2 =(subs(W,{r1x,r1y,r2x,r2y},{r,r,r,r}));
Lorb = imag(q1)*(-epsilon0/k);
S=(k/muu0*w0)*q2;
hw = 2*pi*3*10^8*6.62607015 * 10^-34/(2*pi*lambda);
lorb = hw.*Lorb./S;
plot(r,lorb)

Answer 1

[Walter Roberson]

your values all come out nan. You have terms with really silly exp like

exp(465047415018097955185175015559120439814522795787376405109289013546782463990466768163364960497178088168719688701/146734163107013360162621761657081706047762220035580747699952030138169819136)

which is like exp(3e36)