How do you make a multivariate Newton Raphson method function with N = 50 discrete segments

I have this code here, at the bottom I have a function newrap_me(@hxfunc,J,guessf'). it needs to have N=50 as its discrete segments and perform multivariate Newton Raphson method. How would I go about doing this with the code that I have.

global N;

N = 50;

y=linspace(0,2,N+1);

guess1 = linspace(480,280,N+1);

guessf = [guess1, guess1];

options = optimset('Display','off');

solutions = fsolve(@hxfunc, guessf, options);

hots = solutions(1:N+1);

colds = solutions(N+2:end);

hots = hots';

colds = colds';

fprintf("\n 1")

fprintf("\nPlease Look at Exchanger Length vs. Temperature (fsolve)\n")

figure(1);

plot(y,hots,y,colds);

title('Exchanger Length vs. Temperature (fsolve)');

legend('Hot Stream', 'Cold Stream');

xlabel('Length (m)'); ylabel('T (deg K)');

fprintf("\n 2\n")

J = zeros(2*N+2);

solutions = solutions';

a = 146421;

b = -167.467;

c = 0.721177;

UA = 9600;

for i = 1:N

J(i,i) = a+(b*hots(i))+(c*(hots(i)^2)) - UA;

J(i,i+1) = -1*(a+(b*hots(i+1))+(c*(hots(i+1)^2)));

J(i,N+i+1) = UA;

J(N+i,N+i+1) = (2)*(a+(b*colds(i))+(c*(colds(i)^2)));

J(N+i,N+i+2) = UA - (2*(a+(b*colds(i+1))+(c*(colds(i+1)^2))));

J(N+i,i) = -UA;

end

J(102,102) = 1;

J(101,1) = 1;

fprintf("The norm of the Jacobian is %i", round(norm(J)))

fprintf("\n 3")

ex = newrap_me(@hxfunc,J,guessf');

hots2 = ex(1:N+1);

colds2 = ex(N+2:end);

fprintf("\nPlease Look at Exchanger Length vs. Temperature (NR)\n")

figure(2);

plot(y,hots2,y,colds2);

title('Exchanger Length vs. Temperature (NR)');

legend('Hot Stream', 'Cold Stream');

xlabel('Length (m)'); ylabel('T (deg K)');