MATLAB: Getting error as Matrix dimensions must agree.” can someone help me resolve this

matlab dimensional mismatch error

getting error for my code as matrix dimension must agree "have tried to solve coupled wave equation using ode 45 command " error exaclty is as
Error using .*
Matrix dimensions must agree.
Error in odefun (line 58)
dE_omega_dz(1:length(E_omega)/2) =
fft(1i*conj(E1_t).*E2_t.*exp(1i*Dk.*z)
…
Error in odearguments (line 87)
f0 = feval(ode,t0,y0,args{:}); %
ODE15I sets args{1} to yp0.
Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0,
tfinal, tdir, y0, f0, odeArgs,
odeFcn, …
Error in odefun1 (line 54)
[z, E_omega] =
ode45(@odefun,zspan,E_omega);
code is as- "FUNCTION"
function dE_omega_dz = odefun(z, E_omega,~,~)
dE_omega_dz=zeros(length(E_omega),1);
z=z*10^4;
display(z);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lambda1=800*10^-9;
c=(3*10^8)
lambda2=400*10^-9;
o=1;
y= 28.78076 %22.4431   %22.39002159      %20:0.25:50  22.39002159     22.443 
ne2=1.5687;
no2=1.6934;
no1=1.6614;
r22=2.1*10^-12                   %electro-optic coefficient in m/v
j=1;
L=1
t=1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for Ez=0:500*10^3:85*10^5
noE= no1*(1-0.5*no1^2*r22*Ez)
NEE= ((((sin(y)).^2)/((ne2)^2))+(((cos(y))^2)/((no2)^2)))^-0.5
deltak=-(((4*3.14*(NEE-noE))/ lambda1)) 
Dk(j)=(deltak);
% Dk=deltak;
V=(Ez*4)/10^6
V1(t)=V;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l=800*10^-9;        % lambda
c=30*10^8;
pi=3.1415926535;
omega=2*pi*c/l;
n_2=6.6508*10^-20;
L_NL=6.8286e-18*1.5; 
LGVM=0.6*10^-3;
I0=0.4*10^11;
% you would have to split the fields back in two:
E1_omega = E_omega(1:end/2);
E2_omega = E_omega(end/2+1:end);
% go back to time space to calculate the nonlinear part:
E1_t = ifft(E1_omega);
E2_t = ifft(E2_omega);
figure(786)
x=-300*10^-15:1*10^-15:300*10^-15;
plot(x,fftshift(E1_t.^2));
% pause(0.2)
N=max(size(E1_t));
to=120e-15/1.655; % initial pulse widthin second
dt=1/120e-15;
dw=1/N/dt*2*pi;
%dw=2*pi*c/l;
 w=(-1*N/2:1:N/2-1)*dw;
% and calculate the derivatives:
dE_omega_dz(1:length(E_omega)/2) = fft(1i*conj(E1_t).*E2_t.*exp(1i*Dk.*z) ...
    + 1i*2*pi*n_2*I0*L_NL/l*(abs(E1_t.^2 + ...
    2*abs(E2_t.^2)).*E1_t));
dE_omega_dz(length(E_omega)/2+1:length(E_omega)) = 1i*w.*L_NL/LGVM * E2_t + fft(1i*E1_t.*E1_t.*exp(-1i*Dk.*z)....
    + 1i*4*pi*n_2*I0*L_NL/l*(2*abs(E1_t.^2 + ...
    abs(E2_t.^2)).*E1_t));
j=j+1;
t=t+1;
%k=k+1;
L=L+1;
o=o+1;
end
end
Calling function is as –
%first electric field envelop in fourier space
clc
clear all
close all
Po=0.4*10^11;
% Dk=-60;
Ao=sqrt(Po);
C=0.5;
pi=3.1415926535;
x=-300*10^-15:1*10^-15:300*10^-15;  
c=3e8;
n=1.655;
l=800*10^-9; 
w1=2*pi*c/l;
z=30*10^-3;
beta=-1.5*10^-29;       
i=sqrt(-1);
t=120*10^-15;
Wo=120*10^-15/1.655; %initial pulse width in second
u=Ao*exp(1i*(beta*z-w1*t)).*exp(-((1+i*(-C))/2)*(x./Wo).^2)
 figure(1)
 
plot(x,abs(u).^2,'b');
title('Input Pulse'); xlabel('Time'); ylabel('Amplitude');
hold on
grid on
E11_omega=fft(fftshift(u));
disp(E11_omega)
%second electric field envelop in fourier space
P1=0.4*10^11;
A1=sqrt(P1);
% C=5;
u1=A1*(exp(1i*(beta*z-2*w1*t)).*exp(-((1+i*(-C))/2)*(x./(Wo)).^2));
% figure(2)
% plot(x,abs(u1).^2,'b');
% title('Input Pulse'); xlabel('Time'); ylabel('Amplitude');
% hold on
% grid on
E22_omega=fft(fftshift(u1));
disp(E22_omega)
%  patching of the E1 and E2 electric fields with:
E_omega=[E11_omega,E22_omega]
disp(E_omega)
% figure(3)
% plot(abs(E_omega),'r');
%a=[1:0.1:3]
zspan=[0 30*10^-7];
% Z1=zspan*10;
E_omega_0=[0 0.0625]
[z, E_omega] = ode45(@odefun,zspan,E_omega);
E11_omega = E_omega(:, 1:end/2);
E22_omega = E_omega(:, end/2+1:end);
E1_t = ifft(E11_omega, [], 2);
E2_t = ifft(E22_omega, [], 2);
figure(101)
subplot(2, 2, 1)
[X, Y] = meshgrid(1:size(E11_omega, 2), z);
mesh(X, Y, abs(fftshift(E11_omega, 2)).^2);
subplot(2, 2, 2)
[X, Y] = meshgrid(1:size(E11_omega, 2), z);
mesh(X, Y,abs(fftshift(E22_omega, 2)).^2);
subplot(2, 2, 3)
[X, Y] = meshgrid(x, z);
mesh(X, Y,abs(fftshift(E1_t, 2)).^2);
subplot(2, 2, 4)
[X, Y] = meshgrid(x, z);
mesh(X, Y,abs(fftshift(E2_t, 2)).^2);
% qq(x,:)=E11_omega(end,:);
figure(1)
 plot(x,fftshift(E1_t(end,:)).^2);
% pause(0.1)
% figure(4)
%  plot(z,E_omega,'linewidth',1)
%     xlabel('z')
%    ylabel('phase (radian)')

Best Answer

You have
dE_omega_dz(1:length(E_omega)/2) = fft(1i*conj(E1_t).*E2_t.*exp(1i*Dk.*z) ...
    + 1i*2*pi*n_2*I0*L_NL/l*(abs(E1_t.^2 + ...
    2*abs(E2_t.^2)).*E1_t));
This is in a loop in which you are doing
Dk(j)=(deltak);
where you are incrementing j for each iteration of the for Ez loop. Therefore after the first iteration of for Ez, Dk is not a scalar, and the 601 x 1 expression that you get from 1i*conj(E1_t).*E2_t is being .* by a 1 * j variable. In R2016a and earlier that is an error. In R2016b and later, it is permitted and would mean the same thing as if you had used 
bsxfun(@times, 1i*conj(E1_t).*E2_t, exp(1i*Dk.*z)
producing a 601 x j result -- a size that is increasing in each iteration of the for Ez loop. But the output size dE_omega_dz(1:length(E_omega)/2) is fixed, leading to a problem of trying to output the wrong amount of information.
Perhaps you should be indexing Dk(j) in the statement.
Related Question