Решение систем нелинейных уравнений,метод продолжения,ТЕОРИЯ

function dg=dgxt(x,ps,Eotn)
  for i=1:4;
    xt(i)=x(i);
  end;
  for j=1:4;
    for i=1:4; dx(i)=0; end;
    if abs(xt(j))<1; h=0.1*Eotn;
      else; h=0.1*Eotn*abs(xt(j));
    end;
    dx(j)=h;
    dg(:,j)=(fx(xt+dx)-fx(xt-dx))/(2*h);
  end;
  dg(:,5)=-0.4; dg(5,:)=ps';
end

function p=fp(s,x,Eotn)
  global ps;
  dg=dgxt(x,ps,Eotn);
  b=0;
  b(5)=1;
  p=dg\b';
  norm=sqrt(p'*p); p=p/norm;
  ps=p;
end

function f=fx(x)
  f(1)=2*x(1)+8*x(3)+2*x(4);
  f(2)=2*x(2)+8*x(4);
  f(3)=8*x(1)+2*x(3);
  f(4)=2*x(1)+8*x(2)+2*x(4);
    f(1)=f(1)-6*x(1)^2            +x(1)^3-0.4;
    f(2)=f(2)-6*x(2)^2-4*x(3)*x(4)+x(2)^3-0.4;
    f(3)=f(3)-6*x(3)^2-4*x(2)*x(4)+x(3)^3-0.4;
    f(4)=f(4)-8*x(4)^2-4*x(2)*x(3)+x(4)^3-0.4;
end

fw=fopen('Lr_10.res','w');
  x0=[2.7 5.6 2.3 0.3];
  opt=optimset('TolX',1.0e-5);
  [x,f]=fsolve(@fx,x0,opt);
  fprintf(fw,'\n x='); fprintf(fw,'%10.6f',x);
  fprintf(fw,'\n   f='); fprintf(fw,'%10.6f',f);
fclose(fw);

global ps;
Smax=2500;
Eotn=1.0e-10;
x0=[0; 0; 0; 0; 0];
ps=[0; 0; 0; 0; 1];
opt=odeset('reltol',Eotn);
[s,x]=ode45(@fp,[0 Smax],x0,opt,Eotn);
plot(x(:,1),x(:,5),'-r',x(:,2),x(:,5),'--g',...
            x(:,3),x(:,5),':b',x(:,4),x(:,5),'-.k');
grid on; xlabel('x1, x2, x3, x4'); ylabel('x5 = t');
legend('x1','x2','x3','x4');