title: Best way solving optimization with multiple variables in Matlab?
body: <p>I am trying to compute numerically the solutions for a system of many equations and variables (100+). I tried so far three things: </p>
  
  <ol>
  <li>I now that the vector of p(i) (which contains most of the endogenous variables) is decreasing. Thus I gave simply some starting points, and then was increasing(decreasing) my guess when I saw that the specific p was too low(high). Of course this was always conditional on the other being fixed which is not the case. This should eventually work, but it is neither efficient, nor obvious that I reach a solution in finite time. It worked when reducing the system to 4-6 variables though. </li>
  <li>I could create 100+ loops around each other and use bisection for each loop. This would eventually lead me to the solution, but take ages both to program (as I have no idea how to create n loops around each other without actually having to write the loops - which is also bad as I would like to increase/decrease the amount of variables easily) and to execute. </li>
  <li>I was trying fminsearch, but as expected for that wast amount of variables - no way! </li>
  </ol>
  
  <p>I would appreciate any ideas... Here is the code (this one the fminsearch I tried):</p>
  
  <p><strong>This is the run file</strong>:</p>
  
  <pre><code>clear all
  clc
  
  % parameter
  
  z=1.2;
  w=20;
  lam=0.7;
  tau=1;
  N=1000;
  t_min=1;
  t_max=4;
  M=6;
  a_min=0.6;
  a_max=0.8;
  
  t=zeros(1,N);
  alp=zeros(1,M);
  p=zeros(1,M);
  p_min=2;
  p_max=1;
  
  for i=1:N
  t(i)= t_min + (i-1)*(t_max - t_min)/(N-1);
  end
  
  for i=1:M
  alp(i)= a_min + (i-1)*(a_max - a_min)/(M-1);
  p(i)= p_min + (i-1)*(p_max - p_min)/(M-1);
  end
  
  fun=@(p) david(p ,z,w,lam,tau,N,M,t,alp);
  
  p0=p;
  
  fminsearch(fun,p0)
  </code></pre>
  
  <p><strong>And this is the program-file</strong>:</p>
  
  <pre><code>function crit=david(p, z,w,lam,tau,N,M,t,alp)
  
  X = zeros(M,N);
  pi = zeros(M,N);
  C = zeros(1,N);
  Xa=zeros(1,N);
  Z=zeros(1,M);
  
      rl=0.01;
      rh=1.99;
        EXD=140;
  
        while (abs(EXD)&gt;100)
            r1=rl + 0.5*(rh-rl);  
  for i=1:M
  
  for j=1:N
  
        X(i,j)=min(w*(1+lam), (alp(i) * p(i) / r1)^(1/(1-alp(i))) * t(j)^((z-alp(i))/(1-alp(i))));
        pi(i,j)=p(i) * t(j)^(z-alp(i)) * X(i,j)^(alp(i)) - r1*X(i,j);
  
  end
  
  end
  
          [C,I] = max(pi);
  
          Xa(1)=X(I(1),1);
  
           for j=2:N
      Xa(j)=X(I(j),j);
           end 
  
        EXD=sum(Xa)- N*w;
        if (abs(EXD)&gt;100 &amp;&amp; EXD&gt;0)
            rl=r1;
        elseif (abs(EXD)&gt;100 &amp;&amp; EXD&lt;0)
            rh=r1;
        end
        end
  
  Ya=zeros(M,N);
  
  for j=1:N    
  Ya(I(j),j)=t(j)^(z-alp(I(j))) * X(I(j),j)^(alp(I(j)));
  end
  
    Yi=sum(Ya,2);
  
    if (Yi(1)==0)
            Z(1)=-50;
    end
  
    for j=2:M
        if (Yi(j)==0)
           Z(j)=-50;
        else
        Z(j)=(p(1)/p(j))^tau - Yi(j)/Yi(1);
        end
    end
  
  zz=sum(abs(Z))
  crit=(sum(abs(Z)));
  </code></pre>