Error using mupadmex Error in MuPAD command: The variable is invalid. [stdlib::diff] Error in sym/jacobian (line 34) Jsym = mupadmex('​symobj::ja​cobian',F.​s,x.s);

Dear all My name is Le Thanh Ha. I wrote a code for the monetary research but when I run it, the problem came out.

Please help me find an adequate way to solve this problem

This is my code

      if iBounded==0

      %========== define symbolic variables
      %
      % time index:
      %   0 yesterday
      %   1 today
      %   2 tomorrow
      % define variables used in the model
      syms C0 C1 C2 
      syms Y0 Y1 Y2 
      syms pi0 pi1 pi2 
      syms lambda0 lambda1 lambda2 
      syms No0 No1 No2 
      syms De0 De1 De2 
      syms H0 H1 H2 
      syms R0 R1 R2
      syms s0 s1 s2
      syms w0 w1 w2
      %syms gg0 gg1 gg2
      %syms q0 q1 q2
      %syms x0 x1 x2
      syms pstar0 pstar1 pstar2
      %syms eR0 eR1 eR2                                                           % change this if money rule
      syms em0 em1 em2
      syms gm0 gm1 gm2
      syms m0 m1 m2

      syms logg0 logg1 logg2 
      syms logZ0 logZ1 logZ2 
      syms logTheta0 logTheta1 logTheta2 
      syms logbarPi0 logbarPi1 logbarPi2

      syms gmss0 gmss1 gmss2                                                     % This part for MP rule
      syms pstar_ss0 pstar_ss1 pstar_ss2
      syms sss0 sss1 sss2
      syms wss0 wss1 wss1 
      syms Yss0 Yss1 Yss1 

      syms Clag0 Clag1 Clag2
      syms slag0 slag1 slag2
      %syms rlag0 rlag1 rlag2
      syms mlag0 mlag1 mlag2                                                     % This part for MP rule
      syms gmlah0 gmlag1 gmlag2
      syms logbarPiLag0 logbarPiLag1 logbarPiLag2
      syms pilag0 pilag1 pilag2
      syms Rss0 Rss1 Rss2
      %syms Rsslag0 Rsslag1 Rsslag2
      syms gmss0 gmss1 gmss2
      syms gmsslag0 gmsslag1 gmsslag2

      % define parameters in the model
      bet    = sparam.bet   ;
      gamma  = sparam.gamma ; 
      dn     = sparam.dn    ;
      v      = sparam.v     ; 
      miu    = sparam.miu   ;
      eta    = sparam.eta   ;
      chi    = sparam.chi   ;

      rho_m     = sparam.rho_m;
      s1        = sparam.s1;
      s2        = sparam.s2;

      bar_theta = sparam.bar_theta;
      bar_g  = sparam.bar_g ;                       
      barPiStar = sparam.barPiStar;

      rho_Z     = sparam.rho_Z;
      rho_theta = sparam.rho_theta;
      rho_g     = sparam.rho_g;
      rho_barPi = sparam.rho_barPi;

      sig_Z     = sparam.sig_Z;
      sig_theta = sparam.sig_theta;
      sig_g     = sparam.sig_g ;
      sig_barPi = sparam.sig_barPi;

      Rstar      = @(barPiStar) 1/bet*barPiStar;
      pstar_star = @(barPiStar) ( (1 - eta*barPiStar^((1-chi)*(bar_theta-1)))/(1-eta) )^(1/(1-bar_theta));
      sstar      = @(barPiStar) (1-eta)*pstar_star(barPiStar)^(-bar_theta)/(1 - eta*barPiStar^(bar_theta*(1-chi)));
      wstar      = @(barPiStar) (bar_theta-1)/bar_theta*(1 - eta*bet*barPiStar^(bar_theta*(1-chi)))/(1 - eta*bet*barPiStar^((bar_theta-1)*(1-chi)))*pstar_star(barPiStar);
      No_star    = @(barPiStar) wstar(barPiStar)/(1 - eta*bet*barPiStar^(bar_theta*(1-chi)));
      De_star    = @(barPiStar) 1/(1 - eta*bet*barPiStar^((bar_theta-1)*(1-chi)));
      Ystar      = @(barPiStar) ( (1-gamma*bet)/(dn*(1-gamma))*wstar(barPiStar)/sstar(barPiStar)*bar_g)^(1/(1+v));
      Cstar      = @(barPiStar) 1/bar_g*Ystar(barPiStar);
      lambdastar = @(barPiStar) (1-gamma*bet)/((1-gamma)*Cstar(barPiStar));
      Hstar      = @(barPiStar) Ystar(barPiStar)*sstar(barPiStar);
      mstar      = @(barPiStar) 1/lambdastar(barPiStar)*(Rstar(barPiStar)/(Rstar(barPiStar)-1));
      gmstar     = @(barPiStar) barPiStar;

      %==== I modified so that TFP is stationary (12/20/2012) (zt = 0)
      f1  = lambda1 - 1/(C1 - gamma*Clag1) + gamma*bet/(C2 - gamma*C1);
      f2  = w1 - dn*(H1^v)/lambda1;
      f3  = lambda1 - bet*R1*lambda2/pi2;
      f4  = 1 - (1-eta)*pstar1^(1-bar_theta) - eta*exp(logbarPi1)^(chi*(1-bar_theta)*(1-miu))*pilag1^(miu*chi*(1-bar_theta))*pi1^(bar_theta-1);
      f5  = pstar1 - bar_theta/(bar_theta-1)*No1/De1;
      f6  = No1 - w1 - eta*bet*exp(logbarPi1)^(-(1-miu)*chi*bar_theta)*pi1^(-miu*chi*bar_theta)*pi2^bar_theta*No2;
      f7  = De1 - 1 - eta*bet*exp(logbarPi1)^((1-miu)*chi*(1-bar_theta))*pi1^(miu*chi*(1-bar_theta))*pi2^(bar_theta-1)*De2;
      f8  = s1 - (1-eta)*pstar1^(-bar_theta) - eta*exp(logbarPi1)^(-bar_theta*chi*(1-miu))*pilag1^(-bar_theta*chi*miu)*pi1^(bar_theta)*slag1;
      f9  = Y1 - exp(logZ1)*H1/s1;
      %f10 = R1 - Rss1*(rlag1/Rss0)^rho_R*( (pi1/exp(logbarPi1))^phi_Pi*(Y1/Yss1)^phi_y)^(1-rho_R)*exp(eR1);
      f10 = m1 - (1/lambda1)*(R1/(R1-1));
      f11 = m1 - gm1*mlag1/pi1;
      f12 = gm1 - gmss1*(gmlag1/gmss0)^rho_m*((pi2/exp(logbarPi2))^(-s1)*(Y1/Yss1)^(-s2)*exp(em1))^(1-rho_m);
      f13 = Rss1 - 1/bet*exp(logbarPi1);
      f14 = pstar_ss1 - ( (1 - eta*exp(logbarPi1)^((1-chi)*(bar_theta-1)))/(1-eta) )^(1/(1-bar_theta)); 
      f15 = sss1 - (1-eta)*pstar_ss1^(-bar_theta)/(1 - eta*exp(logbarPi1)^(bar_theta*(1-chi)));
      f16 = wss1 - (bar_theta-1)/bar_theta*(1 - eta*bet*exp(logbarPi1)^(bar_theta*(1-chi)))/(1 - eta*bet*exp(logbarPi1)^((bar_theta-1)*(1-chi)))*pstar_ss1;
      f17 = Yss1 - ( (1-gamma*bet)/(dn*(1-gamma))*wss1/sss1*bar_g)^(1/(1+v));
      f18 = gmss1 - exp(logbarPi1);
      f19 = C1 - 1/exp(logg1)*Y1;
      f20 = Clag2 - C1;
      f21 = gmlag2 - gm1;
      f22 = mlag2 -m1;
      f23 = pilag2 - pi1;
      f24 = slag2 - s1;
      f25 = em2;
      f26 = logg2 - (1-rho_g)*log(bar_g) - rho_g*logg1;
      f27 = logTheta2 - (1-rho_theta)*log(bar_theta)- rho_theta*logTheta1;
      f28 = logZ2 - rho_Z*logZ1;
      f29 = logbarPi2 - (1 - rho_barPi)*log(barPiStar) - rho_barPi*logbarPi1;

      %======================== create function f
      f = [f1;f2;f3;f4;f5;f6;f7;f8;f9;f10;f11;f12;f13;f14;f15;f16;f17;f18;f19;f20;f21;f22;f23;f24;f25;f26;f27;f28;f29];

      % define the vector of controls, y, and states, x
      x  = [Clag1 gmlag1 mlag1 slag1 pilag1 Rss0 wss0 Yss0 pstar_ss0 sss0 em1 logg1 logTheta1 logZ1 logbarPi1];
      y  = [C1 R1 s1 Y1 w1 lambda1 No1 De1 pstar1 pi1 H1];
      xp = [Clag2 gmlag2 mlag2 slag2 pilag2 Rss1 wss1 Yss1 pstar_ss1 sss1 em2 logg2 logTheta2 logZ2 logbarPi2];
      yp = [C2 R2 s2 Y2 w2 lambda2 No2 De2 pstar2 pi2 H2];

      approx=2;
      [fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx]=anal_deriv(f,x,y,xp,yp,approx);

This is my code: "anal_derive.m"

if nargin==5
approx=2;
end

nx = size(x,2);
ny = size(y,2);
nxp = size(xp,2);
nyp = size(yp,2);

n = size(f,1);

%Compute the first and second derivatives of f
fx = jacobian(f,x);
fxp = jacobian(f,xp);
fy = jacobian(f,y);
fyp = jacobian(f,yp);

if approx==2

fypyp = reshape(jacobian(fyp(:),yp),n,nyp,nyp);

fypy = reshape(jacobian(fyp(:),y),n,nyp,ny);

fypxp = reshape(jacobian(fyp(:),xp),n,nyp,nxp);

fypx = reshape(jacobian(fyp(:),x),n,nyp,nx);

fyyp = reshape(jacobian(fy(:),yp),n,ny,nyp);

fyy = reshape(jacobian(fy(:),y),n,ny,ny);

fyxp = reshape(jacobian(fy(:),xp),n,ny,nxp);

fyx = reshape(jacobian(fy(:),x),n,ny,nx);

fxpyp = reshape(jacobian(fxp(:),yp),n,nxp,nyp);

fxpy = reshape(jacobian(fxp(:),y),n,nxp,ny);

fxpxp = reshape(jacobian(fxp(:),xp),n,nxp,nxp);

fxpx = reshape(jacobian(fxp(:),x),n,nxp,nx);

fxyp = reshape(jacobian(fx(:),yp),n,nx,nyp);

fxy = reshape(jacobian(fx(:),y),n,nx,ny);

fxxp = reshape(jacobian(fx(:),xp),n,nx,nxp);

fxx = reshape(jacobian(fx(:),x),n,nx,nx);

else

fypyp=0; fypy=0; fypxp=0; fypx=0; fyyp=0; fyy=0; fyxp=0; fyx=0; fxpyp=0; fxpy=0; fxpxp=0; fxpx=0; fxyp=0; fxy=0; fxxp=0; fxx=0;

end

Calvo_argus line 138 is

[fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx]=anal_deriv(f,x,y,xp,yp,approx);

The problem appear is

Error using mupadmex
Error in MuPAD command: The variable is invalid. [stdlib::diff]

Error in sym/jacobian (line 34)
Jsym = mupadmex('symobj::jacobian',F.s,x.s);

Error in anal_deriv (line 30)
fy = jacobian(f,y);

Error in calvoARsgu (line 138)
[fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx]=anal_deriv(f,x,y,xp,yp,approx);

Error in kkss_ARcalvo (line 44)
[Fconst Pconst Ff Pf Fs Ps] = calvoARsgu(sparam,iBounded);

I greatly appreciate any help.

Thank you