Newton Rhapson for multivariable functions

I am trying to make a program that can solve multivariable equations using the iterative newton rhapson method, but I don't have too much experience in programming or Matlab. This is the formula for the iterative method:

x^(k+1) = x^k - [(f'(x^k))^-1]*f(x^k)

The givens for this are the 2 function f1(x1,x2) and f2(x1,x2) and our initial conditions x1 and x2. Then we have to solve the equations for a certain amount of iterations. x^k, x^(k+1), and f(x^k) are 2x1 matrices. (f'(x^k))^-1 is a 2x2 matrix. Here is what I have so far, but I got stuck with the loop:

%%Newton Rhapson

clear;
clc;
close;

x1 = 0; %initial values (given)
x2 = 1;

initial = (x1,x2);
I = transpose(initial);

f1 = 2x1^2 + x2 - 8; %function 1

f2 = x1^2 - x2^2 + x1*x2 - 4; %function 2

A = (f1,f2)
fx = A.'    % transopose of A to give f(x)

df1 = diff(f1,x1); %differentiates f1 wrt x1)
df2 = diff(f1,x2);
df3 = diff(f2,x1);
df4 = diff(f2,x2);

df = [df1,df2;df3,df4]; %2x2 Matrix with differentials  

for i = 1:4
    y = initial - inv(df)*fx

[Torsten]

Newton iteration is for numerical input, not symbolic.

[curly133]

Right. Edited the original post. Is that better?

[Torsten]

No, "diff" in this case must be applied to a symbolic expression to get the Jacobian.

Try to turn a symbolic "df" into a function handle via "matlabFunction" and continue with pure numerical calculations within the loop.

Best wishes

Torsten.

[curly133]

How do I do this exactly?