初学者请教一个罚函数的问题，我应该怎么改呀？希望各位帮帮忙：下边的这个代码总是运行不出来，会提示出现
Index exceeds matrix dimensions.

Error in Untitled (line 39)
if(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))<=1e-6)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=1e-6)
%罚因子迭代收敛条件
代码如下：
close all
clear all
clc
syms x1 x2 M;              %M为罚因子。
m(1)=1;                    
c=8;                       %c为递增系数。赋初值。
a(1)=0;
b(1)=0;                              
f=x2*(2.56029*x1^2+1638.3447*x1+32190171.4)+38510449*x1^2-2107629.422*x1+9402292899+M*((x1^5-78643.71*x2)^2+(600-x1)^2+(150-x2)^2+(x2-100)^2);%外点罚函数
f0(1)=0;

%求偏导、Hessian元素
fx1=diff(f,'x1');
fx2=diff(f,'x2');
fx1x1=diff(fx1,'x1');
fx1x2=diff(fx1,'x2');
fx2x1=diff(fx2,'x1');
fx2x2=diff(fx2,'x2');

%外点法M迭代循环
for k=1:100                                                         
    x1=a(k);x2=b(k);M=m(k);
    %牛顿法求最优值
    for n=1:100                                                      
        f1=subs(fx1); %求解梯度值和Hessian矩阵
        f2=subs(fx2);
        f11=subs(fx1x1);
        f12=subs(fx1x2);
        f21=subs(fx2x1);
        f22=subs(fx2x2);
        if(double(sqrt(f1^2+f2^2))<=1e-6)                         %最优值收敛条件
            a(k+1)=double(x1);b(k+1)=double(x2);f0(k+1)=double(subs(f));
            break;
        else
            X=[x1 x2]'-inv([f11 f12;f21 f22])*[f1 f2]';
            x1=X(1,1);x2=X(2,1);
        end
    end
if(double(sqrt((a(k+1)-a(k))^2+(b(k+1)-b(k))^2))<=1e-6)&&(double(abs((f0(k+1)-f0(k))/f0(k)))<=1e-6)    %罚因子迭代收敛条件
      %输出最优点坐标，罚因子迭代次数，最优值
      a(k+1)   
      b(k+1)
      k
      f0(k+1)
      break;
else
      m(k+1)=c*m(k);
end
end
% 绘制目标函数曲线图
xx1=0:50:600;
xx2=0:10:150;
for i=1:length(xx1)
    for j=1:length(xx2)
         if( xx1(i)^5-78643.71*xx2(j)<=0&&600-xx1(i)<=0&&150-xx2(j)<=0&&xx2(j)-100<=0 )
            Z(i,j)=xx2(j)*(2.56029*xx1(i)^2+1638.3447*xx1(i)+32190171.4)+38510449*xx1(i)^2-2107629.422*xx1(i)+9402292899;
        else
            Z(i,j)=0;
        end
    end
end
figure(1);
surf(xx2,xx1,Z);
axis([0 600 0 150 0 50000])
title('目标函数曲线图');
xlabel('x1');
ylabel('x2');