complex numbers and Nan's as output using 'fsolve'

 Hi,

I hope someone will spare a few minutes, I am getting complex numbers (NaN's sometimes) as results and wondering if I have incorrect code for the function and/or my use of fsolve is not apporriate here.

clc
clear all

format compact
global a b tau A
I=200;
tau=0.07;
a=0.2;
b=1.5;

Ac=zeros(I+1,1);
Aw=zeros(I+1,1);
Nc=zeros(I,1);
Nw=zeros(I,1);
Theta=zeros(I,1);
Ac0=100;
Aw0=500;
Ac(1,1)=Ac0;
Aw(1,1)=Aw0;
 x0=[0.7,0.3,0.1];
   for i=1:I 
    A=[Ac(i,1),Aw(i,1)];  
    fun=@(x)productivity3(x);
    result=fsolve(fun,x0);
    Nc(i,1)=result(1);
    Nw(i,1)=result(2);
    Theta(i,1)=result(3);
    Ac(i+1,1)=(1+tau*Theta(i,1))*Ac(i,1);
    Aw(i+1,1)=(1+tau*Theta(i,1))*Aw(i,1);

     end
     out=[Nc(I-5:I,1), Nw(I-5:I,1), Ac(I-5:I,1), Aw(I-5:I,1)]

Function

   function F=productivity3(x)      
      global a b tau A

   F=[x(1)+x(2)-1;
       x(1)-(x(3)/a)^(1/b)*(1+tau)*A(1);
       x(2)-(x(3)/a)^(1/b)*(1+tau)*A(2)];

Maybe it is due to choice of parameter values? Thank you. R

Answer 1

[Torsten]

Replace your function by

      function F=productivity3(x) 

      global a b tau A

      F=[x(1)+x(2)-1;

       x(1)-(x(3)^2/a)^(1/b)*(1+tau)*A(1);

       x(2)-(x(3)^2/a)^(1/b)*(1+tau)*A(2)];

This way, it is not possible for fsolve to take the power (^1/b) of a negative number.

Best wishes

Torsten.