I am trying to solve an optimal control problem on the following form of Lorenz system:

function lorenz(du, u, p, t)
theta = p
du[1] = 10.0 * (u[2] - u[1]) + theta[1]
du[2] = u[1] * (28.0 - (u[3])) - u[2] + theta[2]
du[3] = u[1] * u[2] - (8/3) * (u[3]) + theta[3]
end
tspan = (0.0, T) # time span for nonlinear problem
prob = ODEProblem(lorenz,u0,tspan,theta) #ODE problem definition
sol = solve(prob, RK4(),saveat = Δt, dt=Δt) #solving the NL problem

where (theta) is the control parameter. In the previous code, theta is assumed to be constant throughout the time horizon T. However, in an optimal control problem, the control action (theta)
should be updated at each time step. Is there a way to do that without having to manually do the RK4 algorithm?

The differential equations library allows you to define parameters of any type, so functions should also be fine afaik.
Assuming your control parameter only depends on quantities that are known at each step of the integration (i.e time and the state), the following should do the trick:

function controlActionTheta(du,u,t)
Theta = (maximum(u),t,0) #or whatever you want to define there is
return Theta
end
function lorenz(du, u, p, t)
theta = p(du,u,t)
du[1] = 10.0 * (u[2] - u[1]) + theta[1]
du[2] = u[1] * (28.0 - (u[3])) - u[2] + theta[2]
du[3] = u[1] * u[2] - (8/3) * (u[3]) + theta[3]
end
prob = ODEProblem(lorenz,u0,tspan,controlActionTheta)