function [interp,model,s,x,z]=v25_solve(filename,synchron,varnum,usesolver,loop,method,sim1,sim2,sim_start,shock_nodes,approx_order,state_min,state_max,t_convergence,Mu,sigma1,par)
%Consider implementing preallocation!
%!!!Important: Only one of the loops may be used at a time!!!
% There are two cases differentiated by the variable synchron. If
% sychron==1, then two corresponding variables (e.g. the same variable for the two countries) 
% will be varied the same way
% (synchronously) and all the other variables will remain unchanged. If
% synchron==0, then only one variable will be changed and all the rest will
% remain unchanged. If more variables are varied at the same time, the data
% will be stored incorrectly.

%Case 1: Only one variable is varied
if synchron==0;
    for i=1:size(sigma1,3)
     for i2=1:size(Mu,1)
      for i3=1:size(par.ka,1)
       for i4=1:size(par.kb,1)
        for i5=1:size(par.ra,1)
         for i6=1:size(par.rb,1)
          for i7=1:size(par.elastSa,1)
           for i8=1:size(par.elastSb,1)
            for i9=1:size(par.elastDa,1)
             for i10=1:size(par.elastDb,1)
              for i11=1:size(par.gammaa,1)
               for i12=1:size(par.gammab,1)
                for i13=1:size(par.thetaab,1)
                 for i14=1:size(par.thetaba,1)
                  for i15=1:size(par.Pamax,1)
                   for i16=1:size(par.Pbmax,1)
                    for i17=1:size(par.sca,1)
                     for i18=1:size(par.scb,1)
                      for i19=1:size(filename,2)
                        %Calculate a running index
                        sumi=i*i2*i3*i4*i5*i6*i7*i8*i9*i10*i11*i12*i13*i14*i15*i16*i17*i18*i19;
                        totali=size(sigma1,3)*size(Mu,1)*size(par.ka,1)*size(par.kb,1)*size(par.ra,1)*size(par.rb,1)*size(par.elastSa,1)*size(par.elastSb,1)*size(par.elastDa,1)* ...
                            size(par.elastDb,1)*size(par.gammaa,1)*size(par.gammab,1)*size(par.thetaab,1)*size(par.thetaba,1)*size(par.Pamax,1)*size(par.Pbmax,1)*size(par.sca,1)*size(par.scb,1)*size(filename,2);
                        disp(['Solving sub-scenario ' num2str(sumi) ' out of ' num2str(totali) ', i.e. ' num2str(round(sumi/totali*100)) '%'])
                        %Generate RECS model
                        model(1,sumi) = recsmodel(filename(i19).x,struct('Mu',Mu(i2,:),'Sigma',sigma1(:,:,i)^2,'order',shock_nodes),[],struct('eqsolver',usesolver,'display',0));
                        %Replace paramters with new values
                        model(1,sumi).params(1,1)=par.ka(i3);
                        model(1,sumi).params(1,2)=par.kb(i4);
                        model(1,sumi).params(1,3)=par.ra(i5);
                        model(1,sumi).params(1,4)=par.rb(i6);
                        model(1,sumi).params(1,5)=par.elastSa(i7);
                        model(1,sumi).params(1,6)=par.elastSb(i8);
                        model(1,sumi).params(1,7)=1/(1+par.ra(i5));
                        model(1,sumi).params(1,8)=1/(1+par.rb(i6));
                        model(1,sumi).params(1,9)=1/par.elastSa(i7);
                        model(1,sumi).params(1,10)=1/par.elastSb(i8); 
                        model(1,sumi).params(1,11)=par.elastDa(i9);
                        model(1,sumi).params(1,12)=par.elastDb(i10);
                        model(1,sumi).params(1,13)=par.gammaa(i11);
                        model(1,sumi).params(1,14)=par.gammab(i12);
                        model(1,sumi).params(1,15)=par.thetaab(i13);
                        model(1,sumi).params(1,16)=par.thetaba(i14);
                        model(1,sumi).params(1,17)=par.phiamax;
                        model(1,sumi).params(1,18)=par.phibmax;
                        model(1,sumi).params(1,19)=par.Xabmax;
                        model(1,sumi).params(1,20)=par.Xbamax;                    
                        model(1,sumi).params(1,21)=par.Pamax(i15);
                        model(1,sumi).params(1,22)=par.Pbmax(i16);
                        model(1,sumi).params(1,23)=par.sca(i17);
                        model(1,sumi).params(1,24)=par.scb(i18);
                        %Find steady state of model with new parameter
                        %values
                        [model(1,sumi).sss,model(1,sumi).xss,model(1,sumi).zss]=recsSS(model(1,sumi),[],[],struct('eqsolver',usesolver));
                        %Prepare a RECS interpolation structure
                        [interp2,s(:,:,sumi)] = recsinterpinit(approx_order,state_min,state_max);
                        %First Guess
                        [interp2,x(:,:,sumi)] = recsFirstGuess(interp2,model(1,sumi),s(:,:,sumi),model(1,sumi).sss,model(1,sumi).xss,struct('T',t_convergence,'eqsolver',usesolver));
                        %Solve model
                        [interp2,x(:,:,sumi),z(:,:,sumi),~,EXITFLAG] = recsSolveREE(interp2,model(1,sumi),s(:,:,sumi),x(:,:,sumi),struct('eqsolver',usesolver,'loop_over_s',loop,'reemethod',method));
    %OLD:
    %                     if ESIMin==0;
    %                         [SSIM(:,:,:,sumi),XSIM(:,:,:,sumi),ESIM(:,:,:,sumi),STAT(1,sumi),~]=recsSimul(model(1,sumi),interp2,repmat(sim_start,sim1,1),sim2,[],struct('stat',1,'eqsolver',usesolver,'display',0));
    %                     else
    %                         [SSIM(:,:,:,sumi),XSIM(:,:,:,sumi),ESIM(:,:,:,sumi),STAT(1,sumi),~]=recsSimul(model(1,sumi),interp2,repmat(sim_start,sim1,1),sim2,ESIMin(:,:,:,sumi),struct('stat',1,'eqsolver',usesolver,'display',0));
    %                     end;
            %            save(strcat('results/r_',num2str(sigma1(1,1,i)),'_',num2str(sigma1(2,1,i)),'_',num2str(sigma1(2,2,i)),'_',num2str(sigma1(1,2,i)),'.mat'));
            %            moments(:,:,sumi)=STAT.moments;cor(:,:,sumi)=STAT.cor;acor(:,:,sumi)=STAT.acor;
                        interp(1,sumi)=interp2;%interp needs to be assigned this way as it seems that otherwise the recsSolveREE commond does not work properly.
                        if EXITFLAG==0; error('Failure to converge');return;end;
                      end;
                     end;
                    end;
                   end;
                  end;
                 end;
                end;
               end;
              end;
             end;
            end;
           end;
          end;
         end;
        end;
       end;
      end;
     end;
    end;
    
%Case 2: Two corresponding variables will be varied simultaneously 
elseif synchron==1
    for i=1:size(sigma1,3)
     for i2=1:size(Mu,1)
      for i3=1:size(par.ka,1)
       %for i4=1:size(par.kb,1)
        for i5=1:size(par.ra,1)
         %for i6=1:size(par.rb,1)
          for i7=1:size(par.elastSa,1)
           %for i8=1:size(par.elastSb,1)
            for i9=1:size(par.elastDa,1)
             %for i10=1:size(par.elastDb,1)
              for i11=1:size(par.gammaa,1)
               %for i12=1:size(par.gammab,1)
                for i13=1:size(par.thetaab,1)
                 %for i14=1:size(par.thetaba,1)
                  for i15=1:size(par.Pamax,1)
                   %for i16=1:size(par.Pbmax,1)
                    for i17=1:size(par.sca,1)
                     %for i18=1:size(par.scb,1)
                      for i19=1:size(filename,2)
                        %Calculate a running index
                        sumi=i*i2*i3*i5*i7*i9*i11*i13*i15*i17*i19;
                        %Generate RECS model
                        model(1,sumi) = recsmodel(filename(i19).x,struct('Mu',Mu(i2,:),'Sigma',sigma1(:,:,i)^2,'order',shock_nodes),[],struct('eqsolver',usesolver,'display',0));
                        %Replace paramters with new values
                        model(1,sumi).params(1,1)=par.ka(i3);
                        model(1,sumi).params(1,2)=par.kb(i3);
                        model(1,sumi).params(1,3)=par.ra(i5);
                        model(1,sumi).params(1,4)=par.rb(i5);
                        model(1,sumi).params(1,5)=par.elastSa(i7);
                        model(1,sumi).params(1,6)=par.elastSb(i7);
                        model(1,sumi).params(1,7)=1/(1+par.ra(i5));
                        model(1,sumi).params(1,8)=1/(1+par.rb(i5));
                        model(1,sumi).params(1,9)=1/par.elastSa(i7);
                        model(1,sumi).params(1,10)=1/par.elastSb(i7); 
                        model(1,sumi).params(1,11)=par.elastDa(i9);
                        model(1,sumi).params(1,12)=par.elastDb(i9);
                        model(1,sumi).params(1,13)=par.gammaa(i11);
                        model(1,sumi).params(1,14)=par.gammab(i11);
                        model(1,sumi).params(1,15)=par.thetaab(i13);
                        model(1,sumi).params(1,16)=par.thetaba(i13);
                        model(1,sumi).params(1,17)=par.phiamax;
                        model(1,sumi).params(1,18)=par.phibmax;
                        model(1,sumi).params(1,19)=par.Xabmax;
                        model(1,sumi).params(1,20)=par.Xbamax;                    
                        model(1,sumi).params(1,21)=par.Pamax(i15);
                        model(1,sumi).params(1,22)=par.Pbmax(i15);
                        model(1,sumi).params(1,23)=par.sca(i17);
                        model(1,sumi).params(1,24)=par.scb(i17);
                        %Find steady state of model with new parameter
                        [model(1,sumi).sss,model(1,sumi).xss,model(1,sumi).zss]=recsSS(model(1,sumi),[],[],struct('eqsolver',usesolver));
                        %Prepare a RECS interpolation structure
                        [interp2,s(:,:,sumi)] = recsinterpinit(approx_order,state_min,state_max);
                        %First Guess
                        [interp2,x(:,:,sumi)] = recsFirstGuess(interp2,model(1,sumi),s(:,:,sumi),model(1,sumi).sss,model(1,sumi).xss,struct('T',t_convergence,'eqsolver',usesolver));
                        %Solve model
                        [interp2,x(:,:,sumi),z(:,:,sumi),~,EXITFLAG] = recsSolveREE(interp2,model(1,sumi),s(:,:,sumi),x(:,:,sumi),struct('eqsolver',usesolver,'loop_over_s',loop,'reemethod',method));
    %                     if ESIMin==0;
    %                         [SSIM(:,:,:,sumi),XSIM(:,:,:,sumi),ESIM(:,:,:,sumi),STAT(1,sumi),~]=recsSimul(model(1,sumi),interp2,repmat(sim_start,sim1,1),sim2,[],struct('stat',1,'eqsolver',usesolver,'display',0));
    %                     else
    %                         [SSIM(:,:,:,sumi),XSIM(:,:,:,sumi),ESIM(:,:,:,sumi),STAT(1,sumi),~]=recsSimul(model(1,sumi),interp2,repmat(sim_start,sim1,1),sim2,ESIMin(:,:,:,sumi),struct('stat',1,'eqsolver',usesolver,'display',0));
    %                     end;
            %            save(strcat('results/r_',num2str(sigma1(1,1,i)),'_',num2str(sigma1(2,1,i)),'_',num2str(sigma1(2,2,i)),'_',num2str(sigma1(1,2,i)),'.mat'));
            %            moments(:,:,sumi)=STAT.moments;cor(:,:,sumi)=STAT.cor;acor(:,:,sumi)=STAT.acor;
                        interp(1,sumi)=interp2;%interp needs to be assigned this way as it seems that otherwise the recsSolveREE commond does not work properly.
                        if EXITFLAG==0; error('Failure to converge');return;end;
                      end;
                     %end;
                    end;
                   %end;
                  end;
                 %end;
                end;
               %end;
              end;
             %end;
            end;
           %end;
          end;
         %end;
        end;
       %end;
      end;
     end;
    end;        
end;