使用fsolve出现问题，请大家帮我看看，谢谢！

zhdg1989

用fsolve解下边的非线性方程组，调节不同的参数得到的解相同，哪里出错了呢?
function shiyanjidi
clc
syms rho A I E L C1 M K2 l0 C2 m2 l v
rho=7200;
A=4.8;
I=2.5498;
E=5*10^10;
L=40;
C1=0.02;
M=2000;
K2=2392;
l0=6.72;
C2=200;
m2=2390;
l=1;
v=1/2*sqrt(E*I/rho/A)*pi/L;
syms beta f omega omega1 phi gamma1 gamma2 gamma3 mu1 mu2 alpha
beta=A*l^2/8/I
f=2*M*9.8*L^3/(l*E*I*pi^4)
phi=2*M/(rho*A*L)
omega1=sqrt(E*I*pi^4/(rho*A*L^4))
gamma1=C1/(rho*A*omega1)
gamma2=C2/(omega1*m2)
gamma3=K2*rho*A*L^4/(m2*E*I*pi^4)
mu1=2*m2/(rho*A*L)
mu2=1
alpha=l0/l
omega=sqrt(rho*A/E/I)*L*v/pi
p=fsolve('myfun',[0.1 0.1 0.1 0.1 0.1],optimset('Display','off'))
function q=myfun(p)
syms R1 R2;
[K,E]=ellipke(sqrt(p(4)^2/(alpha^2+(p(4))^2)));
R1=2*gamma3*p(4)+8*gamma3*alpha^2*K/pi/p(4)/sqrt(alpha^2+(p(4))^2)-8*gamma3*sqrt(alpha^2+(p(4))^2)*E/pi/p(4);
R2=-2*gamma2*omega*p(4);
q(1)=p(1)+beta*(p(1))^3+3/2*beta*(p(1))*(p(2))^2+3/2*beta*(p(1))*(p(3))^2+omega^2*phi*p(2)-2*f/pi;
q(2)=3*beta*(p(1))^2*p(2)+3/4*beta*(p(2))^3+3/4*beta*p(2)*(p(3))^2+(1-4*omega^2-2*omega^2*phi)*p(2)+2*gamma1*omega*p(3)+mu1*((R1)*cos(p(5))+(R2)*sin(p(5)))+4*f/3/pi;
q(3)=3*beta*(p(1))^2*p(3)+3/4*beta*(p(3))^3+3/4*beta*(p(2))^2*p(3)+(1-4*omega^2-2*omega^2*phi)*p(3)-2*gamma1*omega*p(2)+mu1*((R2)*cos(p(5))-(R1)*sin(p(5)));
q(4)=-4*omega^2*p(4)*cos(p(5))+(R1)*cos(p(5))+(R2)*sin(p(5))+4*omega^2*p(2);
q(5)=4*omega^2*p(4)*sin(p(5))+(R2)*cos(p(5))-(R1)*sin(p(5))+4*omega^2*p(3);
end
end

ChaChing

1.求助完整格式：出错代码和出错提示
2.建议楼主稍微学下函数使用
3.建议数值与符号别混用

ChaChing

1. function shiyanjidi
   
2. %syms rho A I E L C1 M K2 l0 C2 m2 l v
   
3. %syms beta f omega omega1 phi gamma1 gamma2 gamma3 mu1 mu2 alpha
   
4. p=fsolve(@myfun,[0.1 0.1 0.1 0.1 0.1],optimset('Display','off'))
   
5. end
   
6. 
   
7. function q=myfun(p)
   
8. %syms R1 R2;
   
9. rho=7200; A=4.8; I=2.5498; E=5*10^10; L=40; C1=0.02; M=2000; K2=2392; l0=6.72; C2=200; m2=2390; l=1;
   
10. v=1/2*sqrt(E*I/rho/A)*pi/L; beta=A*l^2/8/I; f=2*M*9.8*L^3/(l*E*I*pi^4);
    
11. phi=2*M/(rho*A*L); omega1=sqrt(E*I*pi^4/(rho*A*L^4));
    
12. gamma1=C1/(rho*A*omega1); gamma2=C2/(omega1*m2); gamma3=K2*rho*A*L^4/(m2*E*I*pi^4);
    
13. mu1=2*m2/(rho*A*L); mu2=1; alpha=l0/l; omega=sqrt(rho*A/E/I)*L*v/pi;
    
14. [K,E]=ellipke(sqrt(p(4)^2/(alpha^2+(p(4))^2)));
    
15. R1=2*gamma3*p(4)+8*gamma3*alpha^2*K/pi/p(4)/sqrt(alpha^2+(p(4))^2)-8*gamma3*sqrt(alpha^2+(p(4))^2)*E/pi/p(4);
    
16. R2=-2*gamma2*omega*p(4);
    
17. q(1)=p(1)+beta*(p(1))^3+3/2*beta*(p(1))*(p(2))^2+3/2*beta*(p(1))*(p(3))^2+omega^2*phi*p(2)-2*f/pi;
    
18. q(2)=3*beta*(p(1))^2*p(2)+3/4*beta*(p(2))^3+3/4*beta*p(2)*(p(3))^2+(1-4*omega^2-2*omega^2*phi)*p(2)+2*gamma1*omega*p(3)+mu1*((R1)*cos(p(5))+(R2)*sin(p(5)))+4*f/3/pi;
    
19. q(3)=3*beta*(p(1))^2*p(3)+3/4*beta*(p(3))^3+3/4*beta*(p(2))^2*p(3)+(1-4*omega^2-2*omega^2*phi)*p(3)-2*gamma1*omega*p(2)+mu1*((R2)*cos(p(5))-(R1)*sin(p(5)));
    
20. q(4)=-4*omega^2*p(4)*cos(p(5))+(R1)*cos(p(5))+(R2)*sin(p(5))+4*omega^2*p(2);
    
21. q(5)=4*omega^2*p(4)*sin(p(5))+(R2)*cos(p(5))-(R1)*sin(p(5))+4*omega^2*p(3);
    
22. end

复制代码

p =

    0.0002   -0.0676    0.0467    0.0974   -2.5290