function [k ender]=steepest(f,x,e)

% f=x1-x2+2*x1^2+2*x1*x2+x2^2; 假设f等于
% x=[0;0];
% e=10^(-20);
syms x1 x2 m; %m为学习率
d=-[diff(f,x1);diff(f,x2)] %分别求x1和x2的偏导数，即下降的方向

% - 4*x1 - 2*x2 - 1
% 1 - 2*x2 - 2*x1

flag=1; %循环标志
k=0; %迭代次数

while(flag)
d_temp=subs(d,x1,x(1)) ; %分别对x1，x2求一阶偏导数值
%先把x1带入二阶导函数
d_temp=subs(d_temp,x2,x(2)); %再把x2带入1阶导函数

%以上求出1阶偏导数在x1，x2的函数值

nor=norm(d_temp)%范数 求出此时的偏导数范数

if(nor>=e) %因为在局部最优点偏导数接近0，所以拿来和误差比较

x_temp=x+m*d_temp ; %改变初始点x的值
% 开始迭代x的路径，而且此路径和 该点到数值有关

f_temp=subs(f,x1,x_temp(1)); %将改变后的x1和x2代入目标函数

f_temp=subs(f_temp,x2,x_temp(2))
%分别把更新后的x1，x2带代入 f

h=diff(f_temp,m) %代入含有rate m的x1，x2后的f 对m求导
%此时 是关于m的一元函数，g（m）与f不同，是由（x1，x2）处梯度和f 产生，找到g（m）的最小值，如果最小值不是梯度近似为0，则进行下一次迭代
m_temp=solve(h) ; %求方程，得到当次m 。找出极值
x=x+m_temp*d_temp ;%更新起始点x
%正式求出m并带入 迭代x的方程，求出第二个点
k=k+1 %迭代次数更新
else
flag=0;
end
end
ender=double(x); %此时是满足条件的更新后的 x
end

测试

syms x1 x2;
f=(x1-2)^2+2*(x2-1)^2;  %f可以为任意可微的函数
x=[1;3];
e=10^(-20);
[k ender]=steepest(f,x,e)