Is there a way to export ODE rhs values at solution in addition to solution values

Question

I i'm working on a code that describe biological behavior and i need the data besides the output that the ODE already give, this is better explain below:

function [dy] = diferentialSolver(t,y,cte,VrSTR)
    //[R_pfree]= newtonBissection(cte,y)
    
    [R_pfree] =(cubicRoot(cte,y))
    [cData] = [R_pfree]
    //disp(cubicRoot(cte,y))
    R_pf=[]
    
    // this will select the root in that is a real number, meaning that will always have two complex roots in this system
    if abs(imag(R_pfree(6)))<1D-10 then
        R_pf = abs(real(R_pfree(6)))
        disp("here1")
    elseif (abs(imag(R_pfree(7))))<1.D-10 then
        R_pf = abs(real(R_pfree(7)))
        disp("here2")
    elseif abs(imag(R_pfree(8)))<1.D-10
        R_pf = abs(real(R_pfree(8)))
        disp("here3")
    end

    //disp( newtonBissection(cte,y))
    //redundant to simplify programming
    //[diffdata] = return [R_pfree]
    //R_pf= R_pfree
    V = y(1)
    R_T= y(2)
    // this if statement is to make sure the y(3) can be call either as a matrix or as a single variable.
    if length(y)==3 then
        F = y(3)
    else
        [F]=y(3:$)
    end
    //cte=abs(cte)
    k_fR=cte(1)
    k_fB=cte(2)
    k_dR=cte(3)
    k_dB=cte(4)
    Q = cte(5)
    r_T = cte(6)
    b_T=cte(7)
    K_dQ=cte(8)
    K_dQ2=cte(9)
    K_dR=cte(10)
    K_dB=cte(11)
    r=cte(12) //mi(v)
    k= cte(13)
    vr=VrSTR(:) //promoter strength 
    

    //IGR= 2 Intrinsic Growth Rate, still testing this one.
    //Ext= 1.2    //extinction rate. Not on use yet.
    
    R_gfree =r_T*K_dR/(K_dR+R_pf)
    //disp(R_gfree)    //Free repressor gene concentration at time t
    BM_gfree=b_T*K_dB/(K_dB +R_pf)
    
    //disp(BM_gfree)
   // disp(list(["Rpfree: " + string(R_pfree)],["Rgfree: " + string(R_gfree)], ["BM_gfree: " + string(BM_gfree)]))
    //dy(1)= V*%e^r*(1-V/k)        // Volume over time eq Ricker Model
    
    //the differential equations begins here
    dy(1)= r*V*(1-(V/k))    // Volume over time eq
    
    //does simulate the stabilization of the system well
     dy(2)= (k_fR*vr*R_gfree)-(k_dR*R_pf)-((y(2)/V)*dy(1))    //total repressor over time eq
    // this if statement is to make sure the y(3) can be call either as a matrix or
    //as a single variable.
    if length(y)==3 then
        dy(3)= k_fB*BM_gfree - k_dB*F-F/V*dy(1)    //fluorescence over time
        dy(4)=r_T*K_dR/(K_dR+R_pf)
        dy(5)=b_T*K_dB/(K_dB +R_pf)
    else
        for i = 1:length(F)
            //here we call F as an vector, since we accept that F is a submatrix of y(),
            //containing multiple values of Fluorescence.
            dy(2+i)= k_fB*BM_gfree - k_dB*F(i)-(F(i)/V*dy(1))    //fluorescence over time
            dy(10)=r_T*K_dR/(K_dR+R_pf)
            dy(11)=b_T*K_dB/(K_dB +R_pf)
        end
    end

    disp(dy);
endfunction

I'm trying to receive the dy's 4 and 5, but i only went as far as outputting the last time step, i want all steps, is that even possible? (i mean should be, since disp actualy display the all the dy in time in the console, but scilab apparently can't save as a variable).

I would be in imense appreciation if some one could help me out!

I tried global variables, only got so far as to save the last step in time. I also tried this method above, doesn't even save the other dy's. I tried adding another output in the function like "function [dy,proteinFree] = diferentialSolver(t,y,cte,VrSTR)" didn't work as intended. Also, i tried to save into another function were ODE shouldn't interfere, by calling it inside the diferentialSolver function, but it didn't work. Now i'm out of ideias.

Answer 1

In the actual version of Scilab (2023.1) the fastest way to yield the solution of an ode + the actual rhs values at solution is to reformulate it as a implicit equation with dassl (see https://help.scilab.org/dassl). Here is an example for a simple linear ode:

// original ODE rhs
function yd = rhs(t,y)
  yd = [0 1;-1 0]*y;
end

// equivalent DAE residual
function [r,flag] = res(t,y,yd)
  flag = 0;
  r = yd-rhs(t,y);
end

t = linspace(0,5,100);

// solve with ode function
y0 = [0;1];
y = ode(y0,0,t,rhs);

// solve with dassl and extract y, dy
sol = dassl([y0 rhs(0,y0)],0,t,res);
y = sol(2:3,:);
dy = sol(4:5,:);

You can also call rhs after the ode call, which will involve extra cost if rhs cannot vectorize the computation for multiple y vectors (this is the case for your particular rhs):

y = ode(y0,0,t,rhs);
dy = zeros(y);
for i = 1:size(y,2)
  dy(:,i) = rhs(t(i),y(:,i));
end