下面的matlab程序为何运行不出来？怎么修改？

Fun.m文件
function f=fun( x1,x2 )
f=3*(x1+x2-2)^2+(x1-x2)^2

主程序：
x0=[0.8;0.8];%设置初始点
xk=x0;
ideal_error=10^(-7);
actural_error=1;
d=zeros(2,2);
d(:,1)=[1;0];
d(:,2)=[0;1];
Inc =zeros(2,1);
k=0;
MaxLoopNum=100;%初始化最大迭代次数
while(actural_error>ideal_error&&MaxLoopNum>k)
syms x1;
syms x2;
xktemp =xk;
fun1=fun(x1,x2);
fun1=inline(fun1);
f0=feval(fun1,xk(1),xk(2));%求初始点处函数值
F0=f0;
if k>0
F0=eval(F0);
end
%沿d1方向进行一维搜索
syms a;
syms x1;
syms x2;
xk1=xk+a*d(:,1);
x1=xk1(1);
x2=xk1(2);
fun1=fun(x1,x2);
fxa=diff(fun1,'a');
a=solve(fxa);
xk1=inline(xk1);
xk1=feval(xk1,a);
xk1(1)=eval(xk1(1));
xk1(2)=eval(xk1(2));
syms x1;
syms x2;
fun1=fun(x1,x2);
fun1=inline(fun1);
f1=feval(fun1,xk1(1),xk1(2));
f1=eval(f1);
Inc(1)=f0-f1;
%沿d2方向进行搜索
syms a ;
syms x1;
syms x2;
xk2=xk1+a*d(:,2);
x1=xk2(1);
x2=xk2(2);
fun1=fun(x1,x2);
fxa=diff(fun1,'a');
a=solve(fxa);
xk2=inline(xk2);
xk2=feval(xk2,a);
xk2(1)=eval(xk2(1));
xk2(2)=eval(xk2(2));
syms x1;
syms x2;
fun1=fun(x1,x2);
fun1=inline(fun1);
f2=feval(fun1,xk2(1),xk2(2));
f2=eval(f2);
F2=f2;
Inc(2)=f1-f2;
[Incm,row]=max(Inc);
x3=2*xk2-xk;%计算反射点
syms x1;
syms x2;
fun1=fun(x1,x2);
fun1=inline(fun1);
f3=feval(fun1,x3(1),x3(2));
f3=eval(f3);
F3=f3;
temp1=(F0-2*F2+F3)*(F0-F2-Incm)^2;
temp2=0.5*Incm*(F0-F3)^2; Incm就是Δm
%判断是否更换搜索方向
if(F3<F0&&temp1<temp2)
syms a;
syms x1;
syms x2;
d(:,row)=xk2-xk;
x1=xk(1);
x2=xk(2);
fun1=fun(x1,x2);
fxa=diff(fun1,'a');
a=solve(fxa);
xk=inline(xk);
xk=feval(xk,a);
%不更换搜索方向
else if F2<F3
xk=xk2;
else
xk=x3;
end
end
xkerror=eval(xk2-xktemp);%计算实际收敛精度
actural_error=norm(xkerror);
k=k+1;
end
x=eval(xk);

clc
clear
x0=[0.8;0.8];%设置初始点
xk=x0;
ideal_error=10^(-7);
actural_error=1;
d=zeros(2,2);
d(:,1)=[1;0];
d(:,2)=[0;1];
Inc =zeros(2,1);
k=0;
MaxLoopNum=100;%初始化最大迭代次数
while(actural_error>ideal_error&&MaxLoopNum>k)
syms x1;
syms x2;
xktemp =xk;
fun1=fun(x1,x2);
fun1=inline(fun1);
f0=feval(fun1,xk(1),xk(2));%求初始点处函数值
F0=f0;
if k>0 &&  ~isnumeric(F0)
F0=eval(F0);
end
%沿d1方向进行一维搜索
syms a;
syms x1;
syms x2;
xk1=xk+a*d(:,1);
x1=xk1(1);
x2=xk1(2);
fun1=fun(x1,x2);
fxa=diff(fun1,'a');
a=solve(fxa);
xk1=inline(xk1);
xk1=feval(xk1,a);
xk1(1)=eval(xk1(1));
xk1(2)=eval(xk1(2));
syms x1;
syms x2;
fun1=fun(x1,x2);
fun1=inline(fun1);
f1=feval(fun1,xk1(1),xk1(2));
f1=eval(f1);
Inc(1)=f0-f1;
%沿d2方向进行搜索
syms a ;
syms x1;
syms x2;
xk2=xk1+a*d(:,2);
x1=xk2(1);
x2=xk2(2);
fun1=fun(x1,x2);
fxa=diff(fun1,'a');
a=solve(fxa);
xk2=inline(xk2);
xk2=feval(xk2,a);
xk2(1)=eval(xk2(1));
xk2(2)=eval(xk2(2));
syms x1;
syms x2;
fun1=fun(x1,x2);
fun1=inline(fun1);
f2=feval(fun1,xk2(1),xk2(2));
f2=eval(f2);
F2=f2;
Inc(2)=f1-f2;
[Incm,row]=max(Inc);
x3=2*xk2-xk;%计算反射点
syms x1;
syms x2;
fun1=fun(x1,x2);
fun1=inline(fun1);
f3=feval(fun1,x3(1),x3(2));
f3=eval(f3);
F3=f3;
temp1=(F0-2*F2+F3)*(F0-F2-Incm)^2;
temp2=0.5*Incm*(F0-F3)^2;% Incm就是Δm
%判断是否更换搜索方向
if(F3<F0&&temp1<temp2)
syms a;
syms x1;
syms x2;
d(:,row)=xk2-xk;
x1=xk(1);
x2=xk(2);
fun1=fun(x1,x2);
fxa=diff(fun1,'a');
a=solve(fxa);
xk=inline(xk);
xk=feval(xk,a);
%不更换搜索方向
else if F2<F3
xk=xk2;
else 
xk=x3;
end 
end
xkerror=eval(xk2-xktemp);%计算实际收敛精度
actural_error=norm(xkerror);
k=k+1;
end
x=eval(xk)这是主程序，改正地方如下：

1. % Incm就是Δm前面加%；
2. if k>0 &&  ~isnumeric(F0)
F0=eval(F0);
end加了个~isnumeric(F0)否则如果F0是数值不能用eval
3. 开头加了clc clear可要可不要，建议加

有几段缺了，给你重新写了下，请仔细研读。 clear; clc; load('data'); n=size(x,1); k=size(x,2); x=[ones(n,1) x]; beta_hat=(x'*x)\(x'*y); y_hat=x*beta_hat; e_hat=y-y_hat; d_hat=diag(e_hat.^2); variance_matrix=n*inv(x'*x)*(x'*d_hat*x)*inv(x'*x);