I want to introduce Julia to one course. We have been using MATLAB in this course for many years, which are very robust. For these special problems, I tried Python and it seems not as robust as MATLAB (I am also a beginner of Python), but I managed to solve them all. I heard from friends about Julia, and I tried, but apparently I can hardly solve the same problems (the warning/error is “Instability detected. Aborting”). Here below is the source code, and it seems only the solver `Rosenbrock23`

gives the expected results. This situation makes it difficult for me to explain to others (in MATLAB almost all solvers work, just the efficiency matters).

```
using DifferentialEquations, Plots
function dyth_dtau(dyth, yth, par, t)
# nt = int(length(yth)/2)
nt = par[1] # this can be calculated from the line given above
dz = par[2] # this can be calculated from n (nt)
gamma = par[3]
beta = par[4]
phi = par[5]
Le = par[6]
theta = yth[1:nt]
y = yth[nt+1:2*nt]
rhs = [phi^2 * y[i] * exp(gamma*(1.0 - 1.0/theta[i])) for i in 1:nt]
# dtheta/dtau
dyth[1] = 2 * (theta[2] - theta[1])/dz^2 + beta * rhs[1]
for i in 2:nt-1
dyth[i] = (theta[i+1] - 2 * theta[i] + theta[i-1])/dz^2 + beta * rhs[i]
end
dyth[nt] = (1 - 2 * theta[nt] + theta[nt-1])/dz^2 + beta * rhs[nt]
dyth[1:nt] = dyth[1:nt] / Le
# dy/dtau
dyth[nt+1] = 2 * (y[2] - y[1])/dz^2 - rhs[1]
for i in 2:nt-1
dyth[nt+i] = (y[i+1] - 2 * y[i] + y[i-1])/dz^2 - rhs[i]
end
dyth[2*nt] = (1 - 2 * y[nt] + y[nt-1])/dz^2 - rhs[nt]
end
n = 99
dz = 1/(n+1)
gamma = 30.0
beta = 0.15
phi = 1.1
Le_arr = [1.0, 0.5, 0.25, 0.15, 0.1]
t_span = (0.0, 1.0)
yth0 = ones(2*n)
for i in 1:length(Le_arr)
Le = Le_arr[i]
println(Le)
par = (n, dz, gamma, beta, phi, Le)
prob = ODEProblem(dyth_dtau, yth0, t_span, par)
sol = solve(prob)
if (i == 1)
global plt=plot(sol.t, sol[1, :])
else
global plt=plot!(sol.t, sol[1, :])
end
end
display(plt)
```

Here is the status of the package: " [0c46a032] DifferentialEquations v6.16.0 ".

It is very strange that not all ODE solvers are available. For example: CVODE_BDF() and RadauIIA() are missing.

Could you please give me some advices? Thank you.