fsolve: Index exceeds Matrix dimension

I am trying to solve a system of non-linear equations using fsolve. For that purpose I have defined a function "root2d" of the vector x (which shall represent 15 variables).

The first error I get is " Index exceeds matrix dimensions" with the addition that the error occurs when defining the 15th variable (i.e. vector element) - "Error in root2d (line 19) yd = x(15);". The error is not related to the specific variable but to the 15th position - for whatever reason.

Has anyone had this issue before? Do you have any clue what the reason for this could be?

For your reference there is my code below, the function (system of non-linear equations) as well as the fsolve referring to the function.

function F = root2d(x)

%defining the vector x of 19 endogenous input variables

c       =  x(1);
lambda  =  x(2);
r       =  x(3);
gammma1  =  x(4);
f       =  x(5);
wstar  =  x(6);
Pistarw  =  x(7);
ld      =  x(8);
g1      =  x(9);
mc      =  x(10);
g2      =  x(11);
Pistar  =  x(12);
k       =  x(13);
x       =  x(14);
yd      =  x(15);
vp      =  x(16);
l       =  x(17);
vw      =  x(18);
Rbar  =  x(19);

%defining the parameters

delta   =   0.025;
epsilon =   10;
eta     =   10;
Phi     =   0;
gammma2 =   0.001;
betta   =   0.998;
h       =   0.97;
varpsi  =   8.92;
gammma  =   1.17;
kappa   =   9.51;
alppha  =   0.21;
thetap  =   0.82;
chi     =   0.63;
thetaw  =   0.68;
chiw    =   0.62;
gammmaR =   0.77;
gammmay =   0.19;
gammmaPI=   1.29;
PIbar   =   1.01;
rhod    =   0.12;
rhophi  =   0.93;
sigma_A =   -3.97;
sigma_d =   -1.51;
sigma_phi=  -2.36;
sigma_mu=   -5.43;
sigma_m =   -5.85;
Lambdamu=   3.4e-3;
LambdaA =   2.8e-3;
LambdaYd=   (LambdaA+alppha*Lambdamu)/(1-alppha);
Lambdaz =   LambdaYd;
Lambdac =   LambdaYd;
Lambdaw =   LambdaYd;

%specifying the function

   F1 = (1-h*betta*Lambdaz)/((1-h/Lambdaz)*c)-lambda;
   F2 = Rbar-(PIbar*Lambdaz/betta);
   F3 = r-gammma1;
   F4 = r-(1-(betta/(exp(Lambdaz)*exp(Lambdamu))))*(1-delta)/(betta/(exp(Lambdaz)*exp(Lambdamu)));
   F5 = (1-betta*thetaw*exp(Lambdaz)^(eta-1)*PIbar^(-(1-chiw)*(1-eta)))*f-((eta-1)/eta)*wstar*lambda*PIstarw^(-eta)*ld;
   F6 = (1-betta*thetaw*exp(Lambdaz)^(eta*(1+gammma))*PIbar^(eta*(1-chiw)*(1+gammma)))*f-varpsi*PIstarw^(-eta*(1+gammma))*ld^(1+gammma);
   F7 = (1-betta*thetap*PIbar^((1-chi)*epsilon))*g1-lambda*mc*LambdaYd;
   F8 = (1-betta*thetap*PIbar^(-(1-chi)*(1-epsilon)))*g2-lambda*PIstar*LambdaYd;
   F9 = epsilon*g1-(epsilon-1)*g2;
   F10= k/ld-(alppha/(1-alppha))*(w/r)*exp(Lambdaz)*exp(Lambdamu);
   F11= mc-(1/(1-alppha))^(1-alppha)*(1/alppha)^alppha*w^(1-alppha)*r^alppha;
   F12= (1-thetaw*PIbar^(-(1-chiw)*(1-eta))*exp(Lambdaz)^(-(1-eta)))/(1-thetaw)-PIstarw^(1-eta);
   F13= (1-thetap*PIbar^(-(1-chi)*(1-epsilon)))/(1-thetap)-PIstar^(1-epsilon);
   F14= c+x-yd;
   F15= vp*yd-(exp(LambdaA)/exp(Lambdaz))*k^alppha*ld^(1-alppha)-Phi;
   F16= l-vw*ld;
   F17= (1-thetap*PIbar^((1-chi)*epsilon))*vp/(1-thetap)-PIstar^(-epsilon);
   F18= (1-thetaw*exp(Lambdaz)^eta*PIbar^((1-chiw)*eta))/(1-thetaw)-PIstarw^(-eta);
   F19= k-(exp(Lambdaz)*exp(Lambdamu)*x)/(x*exp(Lambdamu)-(1-delta));

F = [F1; F2; F3; F4; F5; F6; F7; F8; F9; F10; F11; F12; F13; F14; F15; F16; F17; F18; F19];
--------------------------------------------------------------------------------------------------------------
------------------------------------------------SEPERATE CODE-------------------------------------------------
--------------------------------------------------------------------------------------------------------------
%initial guess

c0           = 0.40819;
lambda0      = 2.58777;
r0           = 0.0348988;
gammma10        = 1
f0           = 2.55673;
wstar0       = 1.13352;
PIstarw0     = 1.02002;
ld0          = 0.317845;
g10          = 7.64784;
mc0          = 0.89825;
g20          = 8.4976;
PIstar0      = 1.01869;
k0           = 2.71161;
x0          = 0.0884579;
yd0          = 0.496648;
vp0          = 1.00222;
l0           = 0.318814;
vw0          = 1.00305;
Rbar0           = 0;

x0 = [c0; lambda0; r0; gammma10; f0; wstar0; PIstarw0; ld0; g10; mc0; g20; PIstar0; k0; x0; yd0; vp0; l0; vw0; Rbar0]

fun = @root2d;
F = fsolve(fun,x0);

Answer 1

[Walter Roberson]

You are extracting data from your input variable which is named "x". And then you have

x       =  x(14);

which redefines x to be only whatever was in the 14th input, leaving x a scalar for everything that comes afterwards in the function.

[Daniel Bäumler]

Thanks a lot - it's so obvious but I didn't see it.