Gradient descent residual

Question

I've implemented the gradient descent method for finding roots of a system of nonlinear equations and I am wondering how the residual is determined? Is the residual simply the Euclidean norm (2-norm) which is norm(F,2)? I am trying to plot the residual to see the changes over time but I am quite confused about the term and how you calculate it. If my assumption of norm(F,2) being the residual is correct, I have plotted it and shown it below. If my understand is correct, this shows the algorithm stalling as there is no significant decrease in error each iteration after around iteration 1, resulting in this flat line after the first iteration?

import numpy as np
import array
from numpy.linalg import norm

def gradient_descent(fnon, jac, x0, tol, maxk, maxm, *args):
    k = 0
    x = x0

    F = eval(fnon)(x,*args)

    residuals = []
    arr = []
    arr.append(x)
    
    while (norm(F,2) > tol and k <= maxk):
        J = eval(jac)(x,2,fnon,F,*args)
        delta = -2 * np.matmul(np.transpose(J), F) #
        residuals.append(norm(F,2))
        lamb = 0.01       
        x = x + lamb * delta
        F = eval(fnon)(x,*args)
        arr.append(x)
        k += 1
    if (k >= maxk):
        print('No root found')
    else:
        print('Found root: ')
        print(x)
        
    return(are), residuals

def system(x):
    F = np.zeros((2,1), dtype=np.float64)
    F[0] = x[0]*x[0] + 2*x[1]*x[1] + math.sin(2*x[0])
    F[1] = x[0]*x[0] + math.cos(x[0]+5*x[1]) - 1.2
    return F

def jacob(x,n,fnon,F0,*args):
    J = np.zeros((2,2), dtype=np.float64)
    J[0,0] = 2*(x[0]+math.cos(2*x[0]))
    J[0,1] = 4*x[1]
    J[1,0] = 2*x[0]-math.sin(x[0]+5*x[1])
    J[1,1] = -5*math.sin(x[0]+5*x[1])
    return J

arr, residuals = gradient_descent('system', 'jacob', np.array([[-2],[-1]]), 1e-2, 20, 0)

Plotting norm: