I'm trying to create the Jacobian for a given system of equations f1=s*(y-z), f2=x*(r+z)-y, and f3=x*y-b*z. I am having trouble with my code, I thought this would work but it returns a zeroes matrix.

function[fun]=jacobian(x,y,z,s,r,b)

f1=s*(y-x);
f2=x*(r+z)-y;
f3=x*y-b*z;

delta=.1;

x1=x-delta;
for i=1:3
    f1x(i)=s*(y-x(i));
    f2x(i)=x(i)*(r+z)-y;
    f3x(i)=x(i)*y-b*z;
    x3=x+delta;
end

y1=y-delta;
for i=1:3
    f1y(i)=s*(y(i)-x);
    f2y(i)=x*(r+z)-y(i);
    f3y(i)=x*y(i)-b*z;
    y3=y+delta;
end
z1=z-delta;
for i=1:3
    f1z(i)=s*(y-x);
    f2z(i)=x*(r+z)-y;
    f3z(i)=x*y-b*z(i);
    z3=z+delta;
end

fx1=(f1x(3)-f1x(1))/(2*delta)
fy1=(f1y(3)-f1y(1))/(2*delta);
fz1=(f1z(3)-f1z(1))/(2*delta);
fx2=(f2x(3)-f2x(1))/(2*delta);
fy2=(f2y(3)-f2y(1))/(2*delta);
fz2=(f2z(3)-f2z(1))/(2*delta);
fx3=(f3x(3)-f3x(1))/(2*delta);
fy3=(f3y(3)-f3y(1))/(2*delta);
fz3=(f3z(3)-f3z(1))/(2*delta);

R=[fx1 fy1 fz1;fx2 fy2 fz2;fx3 fy3 fz3]
end

[Walter Roberson]

What are your "x3" and "y3" and "z3" variables for? You overwrite them during each iteration of your loop, but you never use their results for anything.

I notice that in your z loop, your f1z and f2z are independent of indexing with (i), and so are going to be constants ?

[Star Strider]

What elements of {x,y,z,s,r,b} are variables and what are parameters? What are your original equations?