%% N 為積分部數(shù)，h 為積分步長，其值視計(jì)算精度而定

%% x0 為初始值，可任意給定

%% Equfun(x)為非線性方程組表達(dá)式（列向量），

%% Jacobfun(x)為雅可比矩陣

N=100;
h=1/N;
x0=[0 0 0]';
x=x0;

%format long
f=Equfun(x);
b=-h*f;

for i=1:N
A=Jacobfun(x);
k1=inv(A)*b;
A=Jacobfun(x+0.5*k1);
k2=inv(A)*b;
A=Jacobfun(x+0.5*k2);
k3=inv(A)*b;
A=Jacobfun(x+0.5*k3);
k4=inv(A)*b;
x=x+(k1+2*k2+2*k3+k4)/6;
end

disp('The Solution is:')
disp('x=');disp(x);

%%%%%%%%%%%%%%%%%%%%%%%%%

%%%下面是非線性方程組及其雅可比行列式

%%%%%%%%%%%%%%%%%%%%%%%%%

function y=Equfun(x)

%nonliner functions

y=[3*x(1)-cos(x(2)*x(3))-1/2;
x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;
exp(-x(1)*x(2))+20*x(3)+(10*pi-3)/3];


%%%%%%%%%%%%%%%%%%%%%%%%%

function y=Jacobfun(x)

%Jacobi function
%J=zeros(length(x));

y(1,1)=3;
y(1,2)=x(3)*sin(x(2)*x(3));
y(1,3)=x(2)*sin(x(2)*x(3));
y(2,1)=2*x(1);
y(2,2)=-162*(x(2)+0.1);
y(2,3)=cos(x(3));
y(3,1)=-x(2)*exp(-x(1)*x(2));
y(3,2)=-x(1)*exp(-x(1)*x(2));
y(3,3)=20;