%% N 為積分部數，h 為積分步長，其值視計算精度而定
%% x0 為初始值，可任意給定
%% myfun(x)為非線性方程組表達式（列向量），
%% detamyfun(x)為雅可比矩陣
N=100;
h=1/N;
x0=[rand(1) rand(1) rand(1) rand(1)]';
x=x0;
%format long
f=myfun(x);
b=-h*f;
for i=1:N
  A=detamyfun(x);
  k1=inv(A)*b;
  A=detamyfun(x+0.5*k1);
  k2=inv(A)*b;
  A=detamyfun(x+0.5*k2);
  k3=inv(A)*b;
  A=detamyfun(x+0.5*k3);
  k4=inv(A)*b;
  x=x+(k1+2*k2+2*k3+k4)/6;
end
disp('The Solution is:')
format short;
disp('x=');disp(x');
f=myfun(x)'