function f=fxxx(z,x)
global Vz l w0 wy k h H T d Is P I0 p0 M y
l=425.55e-9; %波长；
w0=195e-6; %束腰半径；
k=2*pi/l; %波矢；
h=6.62e-34; %普朗克常量；
H=h/(2*pi);
T=200*10^6*2*pi; %失谐量；
d=5e6*2*pi; %自然线宽；
Is=85; %饱和光强；
%P=2.9e-3; %激光功率；
I0=1.98e5;%8*P/(pi*w0^2); %激光强度；
p0=(I0/Is)*(d^2/(d^2+4*T^2)); %饱和参量常数项；
%p=p0*exp(-2*(y^2+z^2)/w0^2)*sin((2*pi/l*x))^2 %饱和参量；
%U=H*T/2*log(1+p); %光势阱；
M=1.993*10^(-26)*51.991/12; %单个原子的质量/kg；
Vz=926;
E0=M*Vz^2/2;
f(2)=(1+x(2)^2)/(2*(E0-H*T/2*log(1+p0.*exp(-2*(z.^2/w0^2+y.^2/w0^2))*(sin(k*x(1)))^2)))...
*((x(2)*H*T*p0*(sin(k*x(1)))^2.*exp(-2*(z.^2/w0^2+y.^2/w0^2)).*(-4.*z/w0^2))/(2*(1+p0.*exp(-2*(z.^2/w0^2+y.^2/w0^2))*(sin(k*x(1)))^2))...
-H*T*p0.*exp(-2*(z.^2/w0^2+y.^2/w0^2))*sin(2*k*x(1))*k/(2*(1+p0.*exp(-2.*(z.^2/w0^2+y.^2/w0^2))*(sin(k*x(1)))^2)));
f(1)=x(2);
f=f';
return
%球差对沉积的影响。
%close ；
clear;
clc;
global Vz l w0 k h H T d Is I0 p0 M y
l=425.55e-9; %波长；
w0=195e-6; %束腰半径；
k=2*pi/l; %波矢；
h=6.62e-34; %普朗克常量；
H=h/(2*pi);
T=200*10^6*2*pi; %失谐量；
d=5e6*2*pi; %自然线宽；
Is=85; %饱和光强；
%P=2.9e-3; %激光功率；
I0=1.98e5;%8*P/(pi*w0^2); %激光强度；
p0=(I0/Is)*(d^2/(d^2+4*T^2)); %饱和参量常数项；
%p=p0*exp(-2*(y^2+z^2)/w0^2)*sin((2*pi/l*x))^2 %饱和参量；
%U=H*T/2*log(1+p); %光势阱；
M=1.993*10^(-26)*51.991/12; %单个原子的质量/kg；
Vz=926;
E0=M*Vz^2/2;
z0=[-w0,w0];
x0=-l/4:l/200:l/4;
a=p0*k^2*w0^2*H*T/(2*E0);
options=odeset('abstol',1e-30,'reltol',1e-10);
for n=1:length(x0)
[z,x]=ode45('fxxx',z0,[x0(n),0],options);
%计算一个原子在z=0平面上的落点。
%c=length(z);
%re(n)=x(c,1)/0.249e-9;
%ret(n)=x(1,1);
%zp(n)=x0(n)/re(n)
%if abs(re(n))<0.249e-9
%re(n)=0;
%end
%figure(1);
plot(z/w0,x(:,1)/l);
hold on
end
ylabel('纵向位置[w0]','Fontsize',14);
xlabel('入射位置[波长]','Fontsize',14);
编译程序出现错误：
??? Error using ==> mrdivide
Matrix dimensions must agree.
Error in ==> fxxx at 19
f(2)=(1+x(2)^2)/(2*(E0-H*T/2*log(1+p0.*exp(-2*(z.^2/w0^2+y.^2/w0^2))*(sin(k*x(1)))^2)))...
Error in ==> odearguments at 110
f0 = feval(ode,t0,y0,args{:}); %
ODE15I sets args{1} to yp0.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal,
tdir, y0, f0, odeArgs, odeFcn, ...
Error in ==> qiucha at 28
[z,x]=ode45('fxxx',z0,[x0(n),0],options);
------解决方案--------------------------------------------------------
不知道是否已经解决，好像就是个矩阵维数不对应的问题，应该不难！