FEM for nonlinear first-order ODE

Question

Currently I am trying to solve nonlinear Ricatti equation using FEM (Matlab language): $$\frac{d r(z)}{dz} = i k(z) r(z) + \frac{i k(z)}{2}(\epsilon(z) - 1)(1 + r(z))^2$$ $z = [-h/2, h/2]$ and $r(-h/2)=0$. $$k(z) = 2\cdot \pi\cdot n(z)/\lambda$$ $$n(z) = 1 + sin(\pi\cdot(\frac{z}{h} + \frac{1}{2}))$$ $$\varepsilon(z) = n^2(z) $$

First of all I write weak formulation: $$\int\frac{d r(z)}{dz} \cdot N(z) dz= \int i k(z) r(z)\cdot N(z)dz +\int \frac{i k(z)}{2}(\epsilon(z) - 1)(1 + r(z))^2\cdot N(z)dz$$

It should give me system of equations for undetermided coefficients (it is called weighted residuals method as I remember). At this point I have written code for FEM for this problem. It doesnt work correctly, as I am comparing results with ode45. I provide my code below:

clc
clear
h = 10;
n0 = 1;
n_elements = 5;
n_nodes = 2*n_elements+1;
L = h/n_elements;
z = -h/2:L/2:h/2;
tol = 10^(-5);
F = zeros(n_nodes,1);
F(1,1) = 0;
cnt = 0;
n = @(z) 1 + n0*sin(pi*(z/h +0.5));
lambda_range = 0.1*h:0.1*h:10*h;
lda = 10*h;
reflection = zeros(size(lambda_range));
tic
for j = 1:length(lambda_range)
    lda = lambda_range(j);
    while 1
        F1 = assembly(F,n_elements,L,n,lda);
        for i = 1:n_nodes
            if abs(F(i)-F1(i))>tol
                break;
            end
        end
        F = F1;
        cnt = cnt +1;
        if cnt>1000
            break
        end
    end
    reflection(1,j) = (abs(F(end,1)))^2;
end
fprintf('Number of elements=%d\n',n_elements)
fprintf('Number of iterations=%d\n',cnt)

Plotting
plot(lambda_range,reflection,'b','Linewidth',3)
xlabel('\lambda')
ylabel('r(\lambda)')
title('Solution')

plot(z,abs(F).^2,'b','Linewidth',3)
xlabel('z')
ylabel('r(z)')
title('Solution')
toc

%%ODE45 check
tspan = [-0.5 0.5];
y0 = 0;
for j = 1:length(lambda_range)
    lda = lambda_range(j);
    [z,y] = ode45(@(z,y) 1i*2*pi*(1 + sin(pi*(z +0.5)))/lda*y+1i*2*pi*(1 + sin(pi*(z +0.5)))/lda*(((1 + sin(pi*(z +0.5))))^2-1)/2*(1+y)^2, tspan, y0);
    reflection(1,j) = (abs(y(end,1)))^2;
end
plot(z,abs(y).^2,'b','Linewidth',3)
xlabel('z')
ylabel('r(z)')
title('Solution')
toc

plot(lambda_range,reflection,'b','Linewidth',3)
xlabel('z')
ylabel('r(z)')
title('Solution')
toc

%%Functions
function [F1] = assembly(F,n_elements,L,n,lda)
    n_nodes = 2*n_elements+1;
    K_table = [-1/2 -2/3 1/6; 2/3 0 -2/3; -1/6 2/3 1/2];
    K = zeros(n_nodes,n_nodes);
    R = zeros(n_nodes,1);
    lmm = [];
    for i =1:n_elements
        j = (i-1)*2+1;
        lmm = [lmm;[j,j+1,j+2]];
    end


    for i =1:n_elements
        lm = lmm(i,:);

        [int1,int2, int_nl_1,int_nl_2,int_nl_3,f11] = quadraticelement(L,i,n,lda);
        k11 = K_table-1i*(int1+int2) - 1i*(int_nl_1*F(lm(1)) + int_nl_2*F(lm(2)) + int_nl_3*F(lm(3)));
        f1 = 1i*f11;

        K(lm,lm) = K(lm,lm)+k11;
        R(lm)=R(lm)+f1;
    end
    
    %%Boundary condition
    K(1,:) = 0;
    K(1,1) = 1;
    R(1,1) = F(1);

    %%Solution of a system
    d = K\R;
    F1 = d;
end

function [int1,int2, int_nl_1,int_nl_2,int_nl_3,f11, dop, k_table] = quadraticelement(L,i,n,lda)
    gL = [-sqrt(3/5),0,sqrt(3/5)];
    gW = [5/9, 8/9, 5/9];
    k_table = zeros(3); int1 = zeros(3); int2 = zeros(3);int_nl_1 = zeros(3);int_nl_2 = zeros(3);int_nl_3 = zeros(3);f11 = sparse(3,1);
    
    for k = 1:length(gW)
        s = gL(k); w = gW(k);
        N = [(s-1)*s/2, 1-s^2, (s+1)*s/2];
        dns = [(2*s-1)/2, -2*s, (2*s+1)/2];
        coord = [(i-1)*L; (i-1/2)*L; i*L];
        z = L*s/2+(i-1/2)*L;
        J = L/2;
        dz = dns*(1/J);

        f11 = f11 + J*w*N'*(2*pi*n(z)/lda*((n(z))^2 - 1)/2);
        k_table = k_table + J*w*N.'*dz;
        int1 = int1 + J*w*(N')*N*2*pi*n(z)/lda;
        int2 = int2 + J*w*(N')*N*2*pi*n(z)/lda*((n(z))^2 - 1);
        dop = [0, 0, 1];
        
        int_nl_1 = int_nl_1 + J*w*(N')*N*N(1)*(2*pi*n(z)/lda*((n(z))^2 - 1)/2);
        int_nl_2 = int_nl_2 + J*w*(N')*N*N(2)*(2*pi*n(z)/lda*((n(z))^2 - 1)/2);
        int_nl_3 = int_nl_3 + J*w*(N')*N*N(3)*(2*pi*n(z)/lda*((n(z))^2 - 1)/2);
    end
end