Calculate Jacobian for system of nonlinear Equations

Hey guys I'm definitely a beginner in Matlab. And I'm trying to solve the Jacobian for my system of nonlinear equations. I managed to do this symbolically. However, if I try to calculate the Jacobian for certain values, it doesn't work. I always have this Problem

Error using sym/jacobian (line 44) The second argument must be a vector of variables.

However, if I try it for an easier system of equations it works. I don't understand where the problem is. Could please someone help me?

Here is my Code:

if true
 syms x1 x2 x3 x4
%Parameter Values
L=1; Fn=5;Fs=5;bn=0.75;bs=0.75;beta=0.5;rho=0.04;T=0;sigma=1;gamma=1;
x0 = [x1;x2;x3;x4];
%Input Functions
cnn=((Fn*(rho+x2/(1-x1)))/((1-beta)*L*bn)*((sigma-1)/sigma))+sqrt(((Fn*(rho+x2/(1-x1)))/((1-beta)*L*bn)*(sigma-1)/sigma).^2+(Fn*(rho+x2/(1-x1)))/((1-beta)*L*bn)*gamma);
pnn=bn+(Fn*(rho+x2/(1-x1)))/((1-beta)*L*cnn);
%The Functions
f1=x1-(beta*L)/(x2*Fs+bs*L);
f2=(1-beta)*x1*bn*sigma*((pnn/bn).^(1/sigma)*(cnn/sigma+gamma)-gamma)+beta*T-(x3*sigma*((bn/x3).^(1/sigma)*(x4/sigma+gamma)-gamma)*((1-beta)*(1/bn+x1-1)-x2*(Fn/(bn*L))));
f3=(pnn-bn)*(1-beta)*L*cnn-((x3-bn)*L*sigma*((1-beta)*((pnn/x3).^(1/sigma)*(cnn/sigma+gamma)-gamma)+beta*((bn/x3).^(1/sigma)*(x4/sigma+gamma)-gamma)));
f4=(x3-bn)/x3-((((1-beta)*sigma*((pnn/x3).^(1/sigma)*(cnn/sigma+gamma)-gamma)+beta*sigma*((bn/x3).^(1/sigma)*(x4/sigma+gamma)-gamma)))/((1-beta)*(pnn/x3).^(1/sigma)*(cnn/sigma+gamma)+beta*(bn/x3).^(1/sigma)*(x4/sigma+gamma)));
f = [f1;f2;f3;f4];
%initial guess
x0=[0.5,0.5,3,5];
%Jacobian
J=jacobian([f1, f2, f3, f4],[x1, x2, x3, x4])
%calculate Jacobian at initial guess
subs(jacobian([f1, f2, f3, f4],[0.5,0.5,3,5]))

It doesn't look nice, I'm sorry. Hopefully you understand it. I want the Jacobian numerically from the functions f1,f2,f3,f4 at the initial guess. But that doesn't work. I first thought it would be a problem because it is symbolic. However, if I try this simpler version it works. syms x y; jacobian([x - y,x + y.^2]); subs(jacobian([x - y,x + y.^2]),[0,0])

Can someone help me? Thanks very much Barbara

Answer 1

[Torsten]

What about

J=subs(J,[x1,x2,x3,x4],[0.5,0.5,3,5]);

?

Best wishes

Torsten.

Answer 2

[Barbara]

Thank you very much, Torsten, that works!