MATLAB - Newton's Method to solve a system of nonlinear equations

Question

I'm trying to write a function that uses Newton's Method to solve a system of nonlinear equations

Function:

[roots, count, resids, history] = Newtons(func,x0,tol)

with the following inputs and outputs.

I'm struggling with creating the function handle (the first input) and using it in the main function as well as understanding and generating the resids and history outputs. I would really appreciate some help with this function.

What I've done/tried so far: (I'm relatively new to Matlab and believe it is totally wrong and not worth looking at)

% function handle for the first input (func)
function [f, J] = jacobian(f, x)
n=length(x); 
fx=f(x); 
step=1e-10; 

for i=1:n
   xstep = x;
   xstep(i)=x(i)+step;
   J(:,i)=(f(xstep)-fx)/step;
end
end

function [roots, count, resids, history] = Newtons(func,x0,tol)
func = @jacobian;

if nargin == 2
    tol = 10e-10;
end

MAXIT = 50;
xx = x0;
N = 0;
while N < 50
    JJ = func(2);

end