Hi,
If Ex and Ey are the X and Y components and you want to compute the magnitude and phase, you can make some minor changes to your code as follows.
clc
clear all
c = 3*10^8; lambda = 333*10^-9; w1 = c/lambda; w2 = 2*w1;
L=12;M=1000;dx=L/M;x=-L/2:dx:L/2-dx;y=x;
t = 3;
% [X,Y] = meshgrid(x,y);
Ax = cos((pi/180)*(x + y));
Bx = 1i*sin((pi/180)*(x + y));
Ay = cos((pi/180)*(x - y));
By = 1i*sin((pi/180)*(x - y));
Ex = Ax + Bx.*exp(-1i*w1*t);
Ey = Ay + By.*exp(-1i*w2*t);
% Em = [Ex Ey];
figure (1);
magnitude = repmat(abs(Ex).^2,numel(x),1) + repmat((abs(Ey).^2)',1,numel(y));
% imagesc(x,y,abs(Em).^2);
imagesc(x,y,magnitude);
figure (2);
phase = atan(repmat((abs(Ey).^2)',1,numel(y))./repmat(abs(Ex).^2,numel(x),1));
% imagesc(x,y,angle(Em).^2);
imagesc(x,y,phase);