求助，求助，自己的发动机万有特性曲线拟合总出错，望高手指点，谢谢，谢谢！

大树YJS

clear all
%不同转速下的燃油消耗率与扭矩的曲线拟合
T10=[126 102 75 51 26 5];
b1=[322.02 328.27 298.91 317.65 434.86 1500.13];
t1=polyfit(T10,b1,2);
Ttq1=5:120/5:125;
be1=polyval(t1,Ttq1);

T20=[146 121 93 60 30 2];
b2=[310.46 321.76 280.69 285.46 385.25 2670.22];
t2=polyfit(T20,b2,2);
Ttq2=2:145/5:146;
be2=polyval(t2,Ttq2);

T30=[163 134 100 65 32 6];
b3=[295.24 283.98 256.63 283.31 369.67 1204.7];
t3=polyfit(T30,b3,2);
Ttq3=6:157/5:163;
be3=polyval(t3,Ttq3);

T40=[163 129 96 65 32 3];
b4=[285.12 259.86 267.11 288.69 393.4 2310.76];
t4=polyfit(T40,b4,2);
Ttq4=3:160/5:163;
be4=polyval(t4,Ttq4);

T50=[168 134 102 69 32 8];
b5=[276.15 257.5 257.9 297.22 382.34 968.51];
t5=polyfit(T50,b5,2);
Ttq5=8:160/5:168;
be5=polyval(t5,Ttq5);

T60=[177 147 109 68 34 1];
b6=[278.96 265.56 267.22 300.27 397.89 9128.84];
t6=polyfit(T60,b6,2);
Ttq6=1:176/5:177;
be6=polyval(t6,Ttq6);

T70=[176 145 107 72 36 5];
b7=[280.34 266.6 271.75 302.15 405.63 1646.46];
t7=polyfit(T70,b7,2);
Ttq7=5:171/5:176;
be7=polyval(t7,Ttq7);

T80=[176 144 108 72 36 3];
b8=[307.04 276.24 288.11 311.49 411.29 2606.96];
t8=polyfit(T80,b8,2);
Ttq8=3:173/5:176;
be8=polyval(t8,Ttq8);

T90=[167 142 106 72 37 24];
b9=[319.23 282.2 395.9 330.02 423.26 530.78];
t9=polyfit(T90,b9,2);
Ttq9=24:143/5:167;
be9=polyval(t9,Ttq9);

T100=[161 135 105 67 36 18];
b10=[331.79 308.58 324.04 371.34 496.48 747.1];
t10=polyfit(T100,b10,2);
Ttq10=332:415/5:747;
be10=polyval(t10,Ttq10);

T110=[155 129 99 65 38 17];
b11=[355.34 336.81 344.02 386.5 507.92 872.02];
t11=polyfit(T110,b11,2);
Ttq11=355:517/5:872;
be11=polyval(t11,Ttq11);

T120=[145 119 87 59 32 11];
b12=[370.35 364.57 381.07 442.75 623.15 1331.72];
t12=polyfit(T120,b12,2);
Ttq12=11:134/5:145;
be12=polyval(t12,Ttq12);



B=[be1';be2';be3';be4';be5';be6';be7';be8';be9';be10';be11';be12'];
Ttq13=[Ttq1';Ttq2';Ttq3';Ttq4';Ttq5';Ttq6';Ttq7';Ttq8';Ttq9';Ttq10';Ttq11';Ttq12'];
N=[1000*ones(6,1);1500*ones(6,1);2000*ones(6,1);2500*ones(6,1);3000*ones(6,1);3500*ones(6,1);4000*ones(6,1);4500*ones(6,1);5000*ones(6,1);5500*ones(6,1);6000*ones(6,1);6500*ones(6,1)];
G=[ones(72,1),N,Ttq13,N.^2,N.*Ttq13,Ttq13.^2];
A=G\B;               %A为一个6*1矩阵
[n,Ttq]=meshgrid(1000:6500,1:177);% 生成格点矩阵
be=A(1)+n*A(2)+Ttq*A(3)+n.^2*A(4)+n.*Ttq*A(5)+Ttq.^2*A(6);%计算格点上的函数值
P=Ttq.*n/9550;

%外特性输入与样条插值拟合
n1=[983 1501 1998 2500 3003 3497 4000 4500 5003 5502 6010];
T1=[125.245 153.0751 167.7599 169.0451 172.1289 177.784 186.655 188.236 184.9435 178.5386 170.683 158.8832];
n2=1000:6500;  
T2=interp1(n1,T1,n2,'spline');%样条插值拟合并且求解对应n2的转矩T2
TT=repmat(T2,501,1); %TT为501*1401与Ttq同样大小
ii=find(Ttq>TT);     %查找超出外特性转矩的位置
P(ii)=NaN;           %利用非数NaN来裁剪图形
be(ii)=NaN;

%绘制等燃油消耗率曲线和等功率曲线三维拟合图
subplot(1,2,1);
mesh(n,Ttq,P);
hold on;
mesh(n,Ttq,be);
axis([1000,6500,1,177,0,9500]);
xlabel('转速n(r/min)')
ylabel('扭矩Ttq(N.M)')
zlabel('Pe(KW)  be(g/(KW*h))')
title('等燃油消耗率曲线和等功率曲线三维拟合图')

%绘制边界曲线（外特性曲线）
subplot(1,2,2);
plot(n2,T2,'LineWidth',1.5);%生成外特性曲线
axis([1000,6500,1,180]);
xlabel('转速n(r/min)');
ylabel('扭矩Ttq(N*M)');
title('万有特性曲线');
hold on;
%绘制等功率曲线
C=contour(n,Ttq,P,10);   %绘制等功率曲线
clabel(C,'fontsize',8,'color','k'); %标注等位线
hold on;
%绘制等燃油消耗率曲线
C1=contour(n,Ttq,be,15);
clabel(C1,'fontsize',8,'color','k'); %标注等位线
h=legend('外特性曲线','等功率曲线','等油耗曲线');
set(h,'fontsize',8);



??? Error using ==> times
Matrix dimensions must agree.

Error in ==> benteng at 80
G=[ones(72,1),N,Ttq13,N.^2,N.*Ttq13,Ttq13.^2];