矩阵3x3的，又不大，不应该out of memory

请把你的程序打过来,以便诊断.

这是已知矩阵，

A =

[ 418.1625, 12.4299, 0, 0, 0, 0 ]
[ 12.4299, 35.4397, 0, 0, 0, 0 ]
[ 0, 0, 0, 0, 0, 0 ]
[ 0, 0, 0, 12.5800, 0, 0 ]
[ 0, 0, 0, 0, 3.4000*C_55, 0 ]
[ 0, 0, 0, 0, 0, 18.3600]

D =

[ 402.8299, 11.9742, 0, 0, 0, 0 ]
[ 11.9742, 34.1403, 0, 0, 0, 0 ]
[ 0, 0, 0, 0, 0, 0 ]
[ 0, 0, 0, 12.1187, 0, 0 ]
[ 0, 0, 0, 0, 3.2753*C_55, 0 ]
[ 0, 0, 0, 0, 0, 17.6868]

下面是程序，

i=sqrt(-1);
rho=0.0000015;
h=3.4;
f = (1:22); % Frequencies for which experimental results of cg are available
%---------------- -------------------------------------------------------
ka1=sqrt(pi^2/12); % shear correction factor
%--------------------------------------------------------------------------
I0= rho*2*h/2; % intergal of I0 with (-z,z) interval.
I2= rho*2*(h/2)^3/3; % intergal of I2 with (-z,z) interval.
%--------------------------------------------------------------------------
% Symbolic expresion used to perform Algebraic calculation


phi = 0; % Propagation wave angle

for phi_i=phi:phi+1
syms C_55;
syms k ;
syms omega;


c=cos(phi_i/360*2*pi );
s=sin(phi_i/360*2*pi );


if phi_i==0

kx=k*round(c); % wave number in x axis
ky=k*round(s); % wave number in y axis
end

kx=k*(c); % wave number in x axis
ky=k*(s); % wave number in y axis
%---------------------------------------------------------------------------
L33 = -ka1^2*A(5,5)*kx^2-2*ka1^2*A(4,5)*kx*ky-ka1^2*A(4,4)*ky^2+omega.^2.*I0;
L34 = i*ka1*(ka1*A(5,5)*kx+ka1*A(4,5)*ky);
L35 = i*ka1*(ka1*A(4,5)*kx+ka1*A(4,4)*ky);
L44 = D(1,1)*kx^2+2*D(1,6)*kx*ky+D(6,6)*ky^2+ka1^2*A(5,5)-omega.^2.*I2;
L45 = D(1,6)*kx^2+(D(1,2)+D(6,6))*kx*ky+D(2,6)*ky^2+i*ka1^2*A(4,5);
L55 = D(6,6)*kx^2+2*D(2,6)*kx*ky+D(2,2)*ky^2+ka1^2*A(5,5)-omega.^2.*I2;
%--------------------------------------------------------------------------
%Antisymmetric waves matrix:

M = [L33 L34 L35;
L34 L44 L45;
L35 L45 L55];


M=vpa(M);

%-------------------------------------------------------------------------
%Solve the determenent of the matrix and derive the roots

det_M=vpa(det(M));

Sol= solve(det_M,k);

k_A0=sym(Sol(1,1));

omega=2*pi*f;

digits(5);

k_A0=sym((subs(k_A0,omega)));

for p=1:length(f)

Cp_A0((phi_i+1),p)=omega(p)/real(k_A0(p));
end
%--------------------------------------------------------------------------
% Calcualtion of Group velocity values using and analytical differention of
% omega with respect to wavelength

syms omega;

diff_w=diff(Sol(1,1),omega);

omega=2*pi*f;

digits(5);

Cg_A0=subs(diff_w,omega);

end
% --------------------------------------------------------------------------
for m=1:length(f)

Cg_A0_1(:,m)=1/Cg_A0(m);

end

Cg_A0=Cg_A0_1;

end

%Group velocity calculation of Anti-symmetric waves.

Cg2_A0=((Cp_A0(2,:)-Cp_A0(1,:))/(1/360*2*pi));
Cgx_A0=(Cg_A0*cos(phi))-(Cg2_A0.*sin(phi));
Cgy_A0=(Cg_A0*sin(phi))+(Cg2_A0.*cos(phi));
C_55=3.0;
Cgx_A0=subs(Cgx_A0,C_55);
Cgy_A0=subs(Cgy_A0,C_55);
Cg_A0_0=sqrt(Cgx_A0.^2+Cgy_A0.^2);
save Cg_A0_0.mat Cg_A0_0;

我感觉程序没什么问题，但是就是总是超出内存，请高手指点一下，万分感谢！

好久没上网，楼主可能已经解决了吧？下面发表一些自己的看法：
（1）for p=1:length(f)

Cp_A0((phi_i+1),p)=omega(p)/real(k_A0(p));
end
这f的长度本来就是1，何必用for哪？
（2） syms C_55;
这两句必须放到Ａ，Ｄ之前。
（3）vpa求出来的本来就是sym类型的，何必再转换为sym?
希望楼主不要生气，楼主的程序完全没有程序的样子，编程习惯不太好。
楼主可能没学过编程吧，希望您能学一学c语言，学过之后学matlab那是易如反掌；不要以为这是走弯路，这正是一劳永易的
捷径。