FEM for nonlinear first-order ODE

Currently I am trying to solve nonlinear Ricatti equation using FEM:

First of all I write weak formulation:

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)

Number of elements=5

    fprintf('Number of iterations=%d\n',cnt)

Number of iterations=1100

    
    %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

Elapsed time is 1.733690 seconds.

    
    %%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

Elapsed time is 2.092498 seconds.

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

    toc

Elapsed time is 2.272228 seconds.

    
    %%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

Where can be mistake, since for me code is correctly structured and methodology is also correct (or am I wrong?)?

[Torsten]

Don't know about your FEM method, but this should be the result when using ode45.

Up to now, I've not heard about "FEM" in the context of an initial value problem.

h = 10;
zspan = [-h/2 h/2];
r0 = 0;
lambda_range = 0.1*h:0.1*h:10*h;
for j = 1:length(lambda_range)
    lda = lambda_range(j);
    [z,r] = ode45(@(z,r)fun(z,r,h,lda),zspan, r0);
    reflection(1,j) = (abs(r(end,1)))^2;
end
plot(z,abs(r).^2,'b','Linewidth',3)
xlabel('z')
ylabel('r(z)')
title('Solution')

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

function dr = fun(z,r,h,lda)
  n = 1 + sin(pi*(z/h+0.5));
  epsilon = n^2;
  k = 2*pi*n/lda;
  dr = 1i*k*r + 1i*k/2*(epsilon-1)*(1+r)^2;
end

[Andrei]

@Torsten Sorry, but it is not fem part of my code(

[Torsten]

I know, but could you give a link where FEM is used for an initial value problem ?

And I thought you should have at least the reference solution correct ...

I don't want to discourage you, but I don't think somebody here will dive into your FEM code to find a possible mistake.

[Andrei]

I am helpless, like I couldnt find anything similar on github(

FEM is used in assembly function. trial functions are defined inside of quadraticelement function. I have tried to construct matrix of coefficient (it is called K) and rhs R. And since I have nonlinearity I iterate til I obtain answer with exact tolerance.

[Andrei]

It is boundary value problem, and I state zero on a left boundary by setting K(1,:) = 0; K(1,1) = 1; R(1,1) = F(1); After this I find solution F1 = K\R. As I am increasing n_elements results are getting worse and that means that Ive made a mistake inside K matrix. I conducted mathematical analysis and couldnt find mistake

[Torsten]

I am helpless, like I couldnt find anything similar on github(

As I said, I've never heard about using FEM for an initial value problem in one space dimension. Who wants to force you to use this method ?

It is boundary value problem,

No. You have a first-order differential equation with only one condition for r. Such a problem is always an initial value problem - it doesn't matter that the independent variable is named z and not t.

[Andrei]

Oh, I see. I tried to find an example and found it here (not a code and here linear problem is being observed). I tried to apply for my problem and extend it to nonlinear case with simple iterations.

[Andrei]

Who wants to force you,

I study numerical optics this term. Basic task was to obtain solution using Runge-Kutta method of 2nd and 4th order, there was no problems. For "the best" students task was to obtain solution, but with FEM. Since previously I used to work in Fenics package on python (solved scattering from infinite cylinder) I was sure that using FEM will be easy))))))))