function output = Solve_System(N,omega)
    h = 1/N;
    x = [0:h:1]';
    y = [0:h:1]';
    % Solve the nonlinear system by Newton's method
    % The initial value
    f = zeros(N+1);
    u = zeros(N+1);
    for j = 2:N
        for i = 2:N
            f(i,j) = sin(pi*x(i))*sin(pi*y(j));
            u(i,j) = (-1/(2*pi^2))*f(i,j);
        end
    end
    tmax = 20; 
    kmax = 1000; 
    tol = 1.0E-13; % error bound
    F = zeros(N+1);
    t = 0;
    while (t<tmax)
        for j = 2:N
            for i = 2:N
                F(i,j) = (1/h^2)*(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)-4*u(i,j))-u(i,j)^3-f(i,j);
            end
        end
        v = zeros(N+1);
        k = 0;
        rmax = 1;
        % Solve the system by SOR method
        while (k<kmax)&&(rmax>tol)
            rtmp = 0;
            for j = 2:N
                for i = 2:N
                    vGS = ((1/h^2)*(v(i-1,j)+v(i+1,j)+v(i,j-1)+v(i,j+1))+F(i,j))/((4/h^2)+3*u(i,j)^2);
                    vtmp = v(i,j)+omega*(vGS-v(i,j));
                    rtmp = max(rtmp, abs(((4/h^2)+3*u(i,j)^2)*(vtmp-v(i,j))));
                    v(i,j) = vtmp;
                end
            end
            rmax = rtmp;
            k = k+1;
        end
        u = u+v;
        if norm(v)<tol
            break;
        end
        t = t+1;
    end
    % The solution u 
    output = u;
end