To avoid your error message simply specify the derivatives directly (they are easy to obtain from polynomials). With both initial guesses at 1 the system diverges (perhaps that's why you are told only to do 5 iterations!). Setting y = -1 allows the system to converge). Try the following (modify as you see fit):
%%%%% Problem 3.13 %%%%%
xi = 1; % given initial x value
yi = -1; % given initial y value (diverges if both x and y = 1)
i_max = 5; % max number of iterations (use 7 iterations to get f1, f2 = 0)
f1 = @(x,y) -2*x^3 + 3*y^2 + 42; % function 1
f2 = @(x,y) 5*x^2 + 3*y^3 - 69; % function 2
d_f1x = @(x) -6*x.^2; %partial derivative of f1 wrt x
d_f1y = @(y) 6*y; %partial derivative of f1 wrt y
d_f2x = @(x) 10*x; %partial derivative of f2 wrt x
d_f2y = @(y) 9*y.^2; %partial derivative of f2 wrt y
%Jacobian
J = @(x,y) [d_f1x(x) d_f1y(y);
d_f2x(x) d_f2y(y)];
F = @(x,y) [f1(x,y);
f2(x,y)];
for i = 1:i_max
Ji = J(xi, yi);
Fi = F(xi, yi);
del = Ji\Fi;
xy = [xi;yi] - del;
xi = xy(1);
yi = xy(2);
end
disp([xi yi])
disp([f1(xi,yi) f2(xi,yi)])

