用vpasolve求​解非线性方程组时出现​下列问题。

在用vpasolve求解一个非线性方程组时我遇到了以下问题：

程序运行出现以下错误：

错误使用 symengine

Not a square matrix.

出错 sym/privBinaryOp (line 1030)

Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});

出错 ^ (line 330)

B = privBinaryOp(A, p, 'symobj::mpower');

出错 fuzafeixianxingfupin (line 22)

fprintf(fid,'%12.8f %12.8f %12.8f\n',w,sqrt((s.a)^2+(s.b)^2),sqrt((s.c)^2+(s.d)^2));

以下是我的程序：

clc;

clear;

close all;%关闭之前的图形窗口

format long;

global w k ek1 ek2 n f o

k=0.051;ek1=0.00945;ek2=0.0945;n=0.5;f=0.015; o=0.003;

syms a b c d

fid=fopen('fuzafei.data','w');

for w=0:0.01:5

w

eq1=-a.*w.^2+2*w.*ek1.*b+k.*(a-c)+2*ek2.*w.*(b-d)+o.*(3/4*a.^3+3/4*(a.*b.^2))+f;

eq2=-b.*w.^2-2*w.*ek1.*a+k.*(b-d)+2*ek2.*w.*(c-a)+o.*(3/4*b.^3+3/4*(b.*a.^2));

eq3=-c.*n.*w.^2+2*ek2.*w.*(d-b)+k.*(c-a);

eq4=-d.*n.*w.^2+2*ek2.*w.*(a-c)+k.*(d-b);

s=vpasolve([eq1,eq2,eq3,eq4],[a,b,c,d]);

% s.a

% s.b

% s.c

% s.d

fprintf(fid,'%12.8f %12.8f %12.8f\n',w,sqrt((s.a)^2+(s.b)^2),sqrt((s.c)^2+(s.d)^2));

end

fclose(fid);

load fuzafei.data

figure

plot(fuzafei(:,1),fuzafei(:,2),'r-');hold on;

plot(fuzafei(:,1),fuzafei(:,3),'b-');hold on;

请问有没有大神可以帮我看看是哪里出了问题。