Решение системы уравнений в matlab. Не сходится решение

clear
clf
rho=1.225; %Задаём плотность воздуха
%1-номер 2 - контур 1   3 - контур 2 4 - сопротивление1 5 - сопр 2 6 - сопр 3 
%7 - площадь
ksi=[1 3 2 0.5 1.1 0 0.07
    2 1 2 0.7 2.2 1.4 0.1   %x
    3 1 3 1.1 0 0 0.2
    4 2 0 0.15 0.1 0.4 0.09  %k
    5 0 1 0.18 0 0 0.13   %o
    6 4 3 1.13 0 0 0.76
    7 4 5 0.54 0 0 0.77
    8 5 4 0.01 0 0 0.32
    9 3 5 0.66 0 0 0.22
    10 5 3 0.0001 0 0 0.99
    11 4 3 0.22 0 0 0.85
    12 2 4 1 1.1 1.5 1
    13 3 2 1 1.1 5 1
    14 3 0 1 0.1 0.21 1
    15 4 5 1 0 0 1
    16 5 0 1 0 0 1]
 
%1-номер сопр, 2 3 - контуры, 4 напор, 5 расход
HQ=[5 1 0 1000 10;
    3 3 -1 500 6;
    9 5 -3 10000 50;]
%Создаём матрицу расположения напорных элементов
HQ_kont=HQ(:,2:3)
%Создаём матрицу расположения сопротивлений
gg=ksi(:,2:3);
 
schetchik_vetvei=0;
vetv_num=[];
 
%Цикл по созданию матрицы, в которой каждая строка соответствует своей
%ветви. В строке каждой ветви распологаются сопротивления, входящие в неё
for sopr_num=ksi(:,1).'
clear proverka
  
  a_find=ksi(sopr_num,2);
  b_find=ksi(sopr_num,3); 
  
  proverka=find(vetv_num==sopr_num);
  if proverka>0
  1;
  else
      
      vetv_num1=sopr_num;
      schetchik_vetvei=schetchik_vetvei+1;
  for sopr_num2=ksi(sopr_num:max(ksi(:,1)),1).'
      
     if sopr_num2==sopr_num      
        
     elseif ksi(sopr_num2,2)==a_find && ksi(sopr_num2,3)==b_find
vetv_num1=[vetv_num1, sopr_num2.'];
     elseif ksi(sopr_num2,3)==a_find && ksi(sopr_num2,2)==b_find
         
         vetv_num1=[vetv_num1, sopr_num2.'];
     else   
     end 
  end
vetv_num(schetchik_vetvei,1:length(vetv_num1))=vetv_num1;
clear vetv_num1
  end
end
%Сопротивления разложенные по ветвям
disp(vetv_num)
 
%Создание матрицы расположения ветвей между контурами. Значения строк те
%же, что и у vetv_num. Т.е. первая строка vetv_num - сопротивления,
%входящие в условно первую ветвь, а первая строка kont_vetv_num -
%расположение этой ветви между контурами
kont_vetv_num=sort(ksi(vetv_num(:,1),2:3)')'
kont_vetv_num=kont_vetv_num+1;
 
%Использование функции подключённого пакета по теории графов. Создание
%матрицы инцидентности. Количество строк - количество ветвей, по столбцам -
%узлы. Количество узлов на один меньше чем в схеме, но именно столько нужно
%для решения системы Кирхгофа. Последние два столбца матрицы находит
%ошибочно, однако проблема оказывается не в них, да и эти найденные
%столбцы, хоть и не соответствуют узлам, имеют физический смысл
%втекания-вытекания в систему труб. Под комментарием матрица,
%соответствующая действительным узлам, вбитая вручную
Cycles=grCycleBasis(kont_vetv_num)
 
% Cycles =[1     0     0     1     0     0
%         1     1     0     0     0     0
%         1     0     0     0     0     0
%         0     1     0     0     1     0
%         0     1     0     0     0     0
%          0     0     1     1     0     0
%          0     0     1     0     1     0
%          0     0     1     0     0     1
%          0     0     0     1     1     0
%          0     0     0     0     0     1
%          0     0     0     0     1     1];
 
 
%Количество уравнений расходоа=кол-во узлов-1
num_ur_Q=length(Cycles(1,:))
%Количество уравнений напоров
num_ur_H=length(kont_vetv_num(:,1))-num_ur_Q
 
 
%Определение в какой ветви какой напорный элемент установлен
napor_vetv=zeros(1,length(kont_vetv_num'));
for hq=1:length(HQ(:,1))
[napor_str, napor_stolb]=find(vetv_num==HQ(hq,1));
napor_vetv(napor_str)=hq;
end
disp(napor_vetv);
 
%Построение характеристик напорных элементов
Q=[];
H=[];
for i=1:length(HQ(:,1))
Q(i,:)=linspace(0,HQ(i,5),200);
    %Это уравнение характеристики вентилятора. Этот цикл не нужен для
    %программы. Дальше раскроются скобки и в матрицу ZZ добавится множитель
    %Q(i,:)^2, а оставшийся член HQ(i,4) поместится в правую часть
    %уравнения
    H(i,:)=HQ(i,4).*(1-(Q(i,:)./HQ(i,5)).^2);
    figure(1)
    plot(Q(i,:),H(i,:))
    hold on
end
 
%Расчёт Z для каждого сопротивления
Zi=sum(ksi(:,4:6)')'.*rho./2./ksi(:,7).^2
 
% Распределение этих сопротивлений по ветвям
 for vn=1:length(vetv_num(:,1))
Current_Z=0;
for vs=1:length(vetv_num(1,:))
    if vetv_num(vn,vs)~=0
Current_Z=Current_Z+Zi(vetv_num(vn,vs));
    end
end
 
    Z_vetv(vn)=Current_Z;
 
 end 
 disp(Z_vetv)   
 clear Current_Z Zi prov_cords
  
 R=[];
 E=[];
konturi=[]%Это контуры, с которыми работаем по уравнениям напоров
ZZ=zeros(length(kont_vetv_num(:,1)),num_ur_H);
 
%Создание матрицы ZZ для записи уравнений напора в матричном виде
% а также вектора с перечислением контуров. ZZ- аэродинамическое
% сопротивление
for y=2:(num_ur_H+1)
clear R E
konturi=[konturi, y];
 [R,E]=find(kont_vetv_num==y);
 E=ones(size(E)).*(y-1);
for r=1:length(R)
ZZ(R(r),E(r))=Z_vetv(R(r));
end
  end
 
 
 
 
  
 
 
%Создание матрицы HH из 0 и 1 характеризующих мастоположение напорных 
%элементов для удобной записи уравнений напора в матричном виде
HH=zeros(length(HQ(:,1)),num_ur_H);
disp(HQ)
for check_kontur=konturi
    [a,b]=find(abs(HQ(:,2:3))==check_kontur-1);
  if a>0
      for j=1:length(a)
          HH(a(j),check_kontur-1)=(HQ_kont(a(j),b(j)))/(abs(HQ_kont(a(j),b(j))));
      end
  end
end
 
 
disp(ZZ)
%Редактирование матрицы ZZ для учёта действия вентиляторов
for i=1:length(ZZ(1,:))
    for j=1:length(ZZ(:,1))
    if ZZ(j,i)~=0
       for k=1:length(vetv_num(1,:))
           a=find(HQ(:,1)==vetv_num(j,k));
           if a>0
              ZZ(j,i)=ZZ(j,i)+HH(a,i).*HQ(a,4)./HQ(a,5)^2;%Тут прибавляется или вычитается 
              %одно из слагаемых, которое получено раскрытием скобок в уравнении характеристики   
              if ZZ(j,i)<0
                  disp('!!!ВОЗНИКНУТ ОБРАТНЫЕ ТОКИ ЧЕРЕЗ СОПРОТИВЛЕНИЕ №!!!')
                  disp(vetv_num(j,k))
              end
           end
       end
    end
    end
end
disp(ZZ)
 
% 
% 
% for t=1:length(HQ(:,1))
%     for check_kontur=konturi
%         a=find(abs(HQ(t,2:3))==check_kontur-1)>0
%         if a>0
%             
%         if HQ_kont(a)>0
%             HH(t, check_kontur-1)=1;
%         elseif HQ_kont(a)<0
%             HH(t, check_kontur-1)=-1;
%         else
%             disp('ошибка. Рассматривается нулевой контур при создании HH')
%         end
%         end
%     end
% end
% disp(HH)
 
%Выделение матрицы Hmax из исзодной матрицы напорных элементов
 
HHmax=HQ(:,4)'
 
%Этот кусок не нужен
 
    %Правая часть из нулей уравнения напоров в матричной форме
%Answer_H=zeros(1,num_ur_H)
% %   Первое уравнение    
%  Q_find*Cycles = Answer_Q
% 
% %   Второе уравнение
   %  Q_find.^2*ZZ-HHmax*HH
 
%Answer=zeros(1,length(kont_vetv_num(:,1)));
 
%Answer(num_ur_Q+1:num_ur_Q+num_ur_H)=abs(HHmax*HH)
 
%B(1:length(Cycles(:,1)),1:length(Cycles(1,:)))=Cycles;
%B(1:length(Cycles(:,1)),length(Cycles(1,:))+1:length(Cycles(1,:))+length(ZZ(1,:)))=(ZZ)
%HH_v=HHmax*HH;
%x0=Answer/B
 
global num_ur_Q num_ur_H HH_v ZZ Cycles
clear Q
 
fun = @root2d;
x0=1*ones(1,11)
       
E=fsolve(fun,x0)

function Cycles=grCycleBasis(E)
% Function Cycles=grCycleBasis(E) find 
% all independent cycles for a connected simple graph
% without loops and multiple edges
% (fundamental set of circuits).
% For loops and multiple edges you need add new vertexes.
% Input parameter: 
%   E(m,2) - the edges of graph;
%     1st and 2nd elements of each row is numbers of vertexes;
%     m - number of edges.
% Output parameter:
%   Cycles(m,m-n+1) - the Boolean array with numbers of edges.
%     n - number of vertexes;
%     m-n+1 - number of independent cycles.
%     In each column of the array Cycles True value have
%     numbers of edges of this cycle.
% Uses the addition of one edge to the spanning tree and deleting of tails.
% Author: Sergii Iglin
% e-mail: siglin@yandex.ru
% personal page: http://iglin.exponenta.ru
 
nMST=grMinSpanTree(E); % data validation and minimal spanning tree
E=sort(E(:,1:2)')'; % only numbers of vertexes
m=size(E,1); % number of edges
n=max(max(E)); % number of vertexes
Erest=E(setdiff([1:m],nMST),:); % rested edges
nr=m-n+1; % number of rested edges
Cycles=zeros(m,nr); % array for cycles
for k1=1:nr, % we add one independent cycle
  Ecurr=[E(nMST,:);Erest(k1,:)]; % spanning tree + one edge
  A=zeros(n);
  A((Ecurr(:,1)-1)*n+Ecurr(:,2))=1;
  A=A+A'; % the connectivity matrix
  p=sum(A); % the vertexes power
  nv=[1:n]; % numbers of vertexes
  while any(p==1), % we delete all tails
    nc=find(p>1); % rested vertexes
    A=A(nc,nc); % new connectivity matrix
    nv=nv(nc); % rested numbers of vertexes
    p=sum(A); % new powers
  end
  [i1,j1]=find(A);
  incedg=nv(unique(sort([i1 j1]')','rows')); % included edges
  Cycles(:,k1)=ismember(E,incedg,'rows'); % current column
end
return

clear, clc
fun = @root11;
x0 = 10*ones(1,11);  
% x0 = 10*ones(1,3); 
x = fsolve(fun,x0)

function F = root11(Q)
Z1=200;
Z2=100;
Z3=100;
Z4=50;
Z5=100;
Z6=50;
Z7=250;
Z8=100;
Z9=200;
Z10=50;
Z11=100;
 
H5=3000;
Q5=20;
H3=2000;
Q3=40;
H8=4000;
Q8=50;
F(1)=-Q(1)-Q(2)+Q(3);
F(2)=Q(2)-Q(4)+Q(5);
F(3)=-Q(6)-Q(7)+Q(8);
F(4)=Q(1)+Q(6)-Q(9);
F(5)=Q(9)+Q(4)+Q(7)-Q(11);
F(6)=-Q(10)-Q(8)+Q(11);
F(7)=-Z2*Q(2)^2-Z3*Q(3)^2+Z5*Q(5)^2-H5*(1-Q(5)^2/Q5^2)+H3*(1-Q(3)^2/Q3^2);
F(8)=-Z1*Q(1)^2+Z2*Q(2)^2+Z4*Q(4)^2-Z9*Q(9)^2;
F(9)=Z1*Q(1)^2+Z3*Q(3)^2-Z6*Q(6)^2-Z8*Q(8)^2+Z10*Q(10)^2-H3*(1-Q(3)^2/Q3^2)+H8*(1-Q(8)^2/Q8^2);
F(10)=Z6*Q(6)^2-Z7*Q(7)^2+Z9*Q(9)^2;
F(11)=Z7*Q(7)^2+Z8*Q(8)^2+Z11*Q(11)^2-H8*(1-Q(8)^2/Q8^2);
 
% F(1)=Q(1)-Q(2)-Q(3);
% F(2)=-200*Q(2)^2+350*Q(3)^2+1000*(1-Q(2)^2/10^2)-500*(1-Q(3)^2/10^2);
% F(3)=170*Q(1)^2+200*Q(2)^2-500*(1-Q(1)^2/7^2)-800*(1-Q(3)^2/10^2);
end