function Newton()
x0=[0.1;0.5];
x1=x0-inv(myJacobi(x0))*myfun(x0);
while norm(x1-x0)>1e-3
x0=x1;
x1=x0-inv(myJacobi(x0))*myfun(x0);
end
x1
msgbox('MATLAB编程答疑，请加QQ: 1530497909','MATLAB答疑','help')
web http://url.cn/NSFcAs -browser
end
function f=myfun(x)
syms x1 x2
f1=(15*x1+10*x2)-((40-30*x1-10*x2)^2*(15-15*x1))*5e-4;
f2=(15*x1+10*x2)-((40-30*x1-10*x2)*(10-10*x2))*4e-2;
f=[f1;f2];
x1=x(1);
x2=x(2);
f=eval(f);
end
function J=myJacobi(x)
syms x1 x2
f1=(15*x1+10*x2)-((40-30*x1-10*x2)^2*(15-15*x1))*5e-4;
f2=(15*x1+10*x2)-((40-30*x1-10*x2)*(10-10*x2))*4e-2;
J=[diff(f1,x1) diff(f1,x2)
diff(f2,x1) diff(f2,x2)];
x1=x(1);
x2=x(2);
J=eval(J);
end
