Implementation of Integral Cost function in Matlab

Please, I need ideas how to simulate the model in the attached document. writing code for equation 39 - 41 is trivial, however, I am not sure how to write the code for equation 42.

I implemented the code below for one time step and assuming that the control input u = 0.1. The question is it correct to compute the optimal cost function like this or there is a better way. Please, find attached the model.

tspan = [0,1];
x0 = 1;
u = 0.1;
tau = 1;
xn = 1;

[time, dxdt, J] = plant_dynamics(tspan,x0,u, tau, xn);

%% System Dynamics
function [time, dxdt, J] = plant_dynamics(tspan,x0,u, tau, xn)

    [time, dxdt] = ode23(@solve_ode,tspan,x0);
    
    J = xn + integral_cost(dxdt, u);

    function dx = solve_ode(t,x)

        A = 1 + tau/12000;
        B = 1 + 0.25 * sin(2*pi*t/3000); 
        dx = A*x + B * u;
    
    end

    function int_J = integral_cost(dxdt, u)
        x = dxdt;
        
        integral_J = x.^2 + u.^2; 
        
        int_J = trapz(integral_J); 
        
    end

end

Answer 1

[VBBV]

tspan = [0,1];
x0 = 1;
u = 0.1;
tau = 1;
xn = 1; % noise 
[time, dxdt, J] = plant_dynamics(tspan,x0,u, tau, xn);
subplot(211)
plot(time,dxdt); title('Plant response')
subplot(212)

plot(time,J);title(' cost function (J)  varying with noise input ')

%% System Dynamics
function [time, dxdt, J] = plant_dynamics(tspan,x0,u, tau, xn)

    [time, dxdt] = ode23(@solve_ode,tspan,x0);
    
    for k = 1:length(dxdt)
        
        J(k,:) = (rand(1)*xn *dxdt).^2 + integral_cost(dxdt, u,xn);  %  add noise here 
    end 

    function dxdt = solve_ode(t,x)

        A = 1 + tau/12000;
        B = 1 + 0.25 * sin(2*pi*t/3000); 
        dxdt = A*x + B * u  ;
    
    end

    function int_J = integral_cost(dxdt, u,xn)
        
        x = dxdt;
        
        integral_J = x.^2 + u.^2; 
        
        int_J = trapz(integral_J,xn); 
        
    end

end