I am trying to solve a ODE with ModelingToolkit.jl, unfortunately I have been waiting for the result for more than 11 hours.
I would to use Multithread, or some CUDA package in order to increase this awful performance, but I am pretty new to this (not only to Julia).

using ModelingToolkit, DifferentialEquations;
using Plots;
using DataFrames, CSV;

@parameters t ξ; @variables φ(t) ψ(t);

N = 2;
n = 4;
V₀ = 0.002;
B = φ^N;
V = V₀ * φ^n;

D = Differential(t);

dB = expand_derivatives(D(B)/D(φ));
dV = expand_derivatives(D(V)/D(φ));
ddB = expand_derivatives(D(dB)/D(φ));
ddV = expand_derivatives(D(dV)/D(φ));

There is a few reasons why your code is running slow. Mainly that all your variables are in the global scope and so the compiler can not generate efficient code for it. The first thing I’d do is wrap your entire code into a function. That should help a lot. Then consider reading the Performance Tips.

That being said, I am not sure why that piece of code is taking more than 11 hours to run. What part of the code is getting hung?

Lastly, please provide julia version versioninfo() and package versions by julia> ] st

I do not think it is because of that. The ModelingToolkit converts system of equations into a function by itself and calls with given initial state u01 and p1. @leonandonayre you can try a small time step(if it is oscillating fast then use a fraction of the period of oscillation) for tspan to see if it really shows something or if it stuck in somewhere in compilation. If it still takes much time then interrupt and show us some stacktrace from the interruption.

The problem appears to be that the model completely blows up between x=14 and x=15. If you use Rodas5() as your solver, you will see that it finishes in about 5 seconds, but stops after around x=14 because the solution blows up to infinity.

You see that every robust solver from each language all blows up to infinity at almost precisely the same point in time. That is very strong evidence that the true solution of the model you had written down goes to infinity at that time.

It seems like the automatic stiffness detection falls into a loop (which we should fix), but at least your issue is that the differential equation you wrote down is divergent and you should probably double check for sign issues.