Problem simulating Matlab ODE15s in Julia

Here is the code I have,

using Pkg
using PyCall


using CSV; using DataFrames; using Dierckx;using DifferentialEquations;using  LSODA;using MATLABDiffEq; using ParameterizedFunctions;using IfElse;

p = [
	1.35e-2  #k1
    6.33e-1  #k2
    5.00e-5  #k3
    1.00e-3  #k4
    3.80e-3  #k5
    5.82e-1  #k6
    2.20e-2  #k7
    4.71     #k8
    1.08e-2  #k9
	#2.6     REMOVED
    1.35     #sigma /index 10
    0.63     #Km /index 11
    5.0      #Gb /index 12
    8.58     #Ib /index 13
]    

input = [
        75000, # D
        70, # Mb
        0 # t_meal_start
    ]

c = [0.00551, 
    17/70, 
    0.043,
    1,
    31, 
    3, 
    13/70, 
    9, 
    30, 
    0.1, 
    0.043*(p[11]+p[12])/p[12]-p[5]*1*p[13], 
    p[7]*p[12]/(1*31*p[13])*30 
] 
function differentialequations_1(du,u,p,t)
    Mg,Gpl,Ipl,Int=u
   
    mg = IfElse.ifelse(t> input[3], p[10].*p[1].^p[10].*(t-input[3]).^(p[10]-1) .* exp(-(p[1].*(t-input[3])).^p[10]) .* input[1],0)
    
    mgpl = p[2] .* Mg
    du[1] = mg - mgpl
    
   
    gli    = c[3] - p[3] .* (Gpl-p[12])-p[4] .* c[4] .* (Ipl-p[13])

    ggu    = p[2] .* (c[1]/(c[2]*input[2])) .* Mg;
    gnoni  = c[11] .* (Gpl./(p[11]+Gpl));
    gi     = p[5] .* c[4] .* Ipl .* (Gpl./(p[11]+Gpl))
    gre= IfElse.ifelse(Gpl > c[8], c[10] ./ (c[2]*input[2]) .* (Gpl - c[8]), 0 )
    du[2] = gli + ggu - gnoni - gi - gre

    
    
    t_lowerbound= t- c[9];

    if (time > c[9]) && (length(t_saved)>1) && (length(t_saved) == length(Gpl_saved))
        spl=Spline1D(t_saved,Gpl_saved)
        Gpl_lowerbound= spl(t_lowerbound);
    else
        Gpl_lowerbound = Gpl_saved[1]; 
    end 

    du[4] = (Gpl-p[12]) - (Gpl_lowerbound-p[12]);
    ipnc = (c[4].^-1) .* (p[6].*(Gpl-p[12]) +(p[7]/c[5]).*Int + (p[7]/c[5]) .*p[12]* c[9] + (p[8].*c[6]) .* du[2])

    iliv    = c[12].*Ipl;
    iif     = p[9].*(Ipl-p[13])
    
    du[3] = ipnc - iliv - iif

end



function compute(u0, local_G_saved,local_I_saved)
    
    global time=0
    global t_saved=[0]
    
    convert(Any, local_G_saved)
    convert(Any, local_I_saved)
    while   time < 120
        tspan = (time, time+1)
        
        prob = ODEProblem(differentialequations_1,u0,tspan,p)
        sol=solve(prob, MATLABDiffEq.ode15s())
        time=  time+1
        convert(Any, sol)
        push!(t_saved, time)
        push!(local_G_saved, last(sol[2,:]))
        push!(local_I_saved,last(sol[3,:]))
        u0=last(sol)
        
        
    end
    
    return local_G_saved, local_I_saved, u0
end 

global Gpl_saved= 5.0
global Ipl_saved=11.5
u0=[0.0,5.0,11.5,0.0]
p[1]=0.00933
p[5]=0.01
p[6]=0.3
p[8]=3.3
simulated_G_saved, simulated_I_saved, u= compute(u0, Gpl_saved,Ipl_saved);

which when run provides the error

Unable to resolve the name mxsol.stats.

Now, I have tried reproducing the error is a smaller version but unfortunately I could not. With every attempt there were some other issues popping up. In some cases, the 4th u variable entirely disappeared in the 2nd iteration. Please run the code and see if you can reproduce the error and see where the flaw is.
I tried to follow the julia recommended coding practices, but please feel free to provide additional pointers if necessary.

Thanks in advance :slight_smile:

Better:

using CSV, DataFrames, Dierckx, DifferentialEquations, LSODA, MATLABDiffEq, ParameterizedFunctions, IfElse

OK, I still have to look at the main issue…

Just realized that you want to use Matlab to solve a differential equation. I do not have a license, so I cannot help you, sorry… Why do you want to use Matlab in the first place?

I noticed a difference between the simulated results from Matlab and Julia. So by using the MatlabDiffEq package I want to ensure my Julia implementation is actually correct.

What version of MATLAB is that?

Matlab R2021b academic use

Try this branch: Add safety around the solver stats by ChrisRackauckas · Pull Request #31 · SciML/MATLABDiffEq.jl · GitHub

Thanks for what seems to be a patch, but could someone please tell me how I can try out this patch. Can I just modify MATALABDiffEq.jl locally or do I have to pull the entire branch? Or is there a different procedure?

On the phone now, but i think pkg> add MATLABDiffEq#ChrisRackauckas-patch-1 would add that specifik version.