I’m working in a mathematical model developed using Julia. I’ve the same model done in Mathematica and, strangely, the solution for the system of ODEs takes more time in Julia than in Mathematica (around 60s vs 24s).

As it’s my first model in Julia I think I might have built the model using some non-efficient ways of implementation.

Here is a sample of my code:

```
#Some data from spreadsheet
tabs=(Dates.value.(data[:,1])-Dates.value(data[1,1])*ones(size(data)[1]))*24*60*60; #Time (s)
P=Interpolations.LinearInterpolation(tabs,data[:,13]*0.001/(24*60*60)); #Rainfall (m/s)
Qin=Interpolations.LinearInterpolation(tabs,data[:,2]*1/(24*60*60)); #Influent water flow rate (m3/s)
ET0=Interpolations.LinearInterpolation(tabs,data[:,14]*0.001/(24*60*60)); #Reference evap019 otranspiration (m/s)
T=Interpolations.LinearInterpolation(tabs,data[:,3]+273.15*ones(size(data)[1])); #Average temperature (K)
fp=Interpolations.LinearInterpolation(tabs,data[:,12]/24); #Photoperiod
I0=Interpolations.LinearInterpolation(tabs,data[:,9]*1000000/(24*60*60)); #Average incident daylight intensity (W/m2)
w=Interpolations.LinearInterpolation(tabs,float.(data[:,10])); #Average wind speed (m/s)
wmax=Interpolations.LinearInterpolation(tabs,float.(data[:,11])); #Maximum wind speed (m/s)
#Some numerical parameters
g=9.8; #Standard gravity of Earth (N/kg)
hs=0.45; #Piezometric head at the outlet (m)
d=(200-2*4.9)*10^(-3); #Outlet's pipe diameter (m)
e=0.0015*10^(-3); #Pipe's absolute roughness (m)
f=0.015; #Outlet pipe's friction factor
L=6.00; #Outlet's pipe lenght (m)
KL=0.75; #Outlet's pipe kinetic energy factor due to minor losses at the elbow
n=3; #Number of pipes at the outlet
S0=100000; #Reference wetland area (m2)
h0=0.50; #Reference wetland depth (m)
p=4.55; #Shape parameter
Kg=10.0*0.001/(24*60*60); #Soil's hydraulic conductivity constant (m/s)
Kc=1.2; #Vegetation coefficient
F=80; #Fetch (m)
RF=1:10; #Resuspension factors for different groups of birds
Ne=[15,50,0,300,0,40+50,40,0,0,0]; #Number of birds of each type
Xin=[1,0,3.1,3.1,12,12,0.007,0.001,0.001,12,12,0,0,18.3,1,4.5]*1000/1000000; #Input concentrations (kg/m3)
X0=[1,1800,3.1,3.1,12,12,0.007,0.001,0.001,12,12,4505.4,0.005,18.3,1,4.5]*1000/1000000; #Initial conditions
#Some functions
w0(t)=117.1267*h(t)/(pi*tanh(pi))*(g/F)^0.5; #Wind speed limit for resuspension (m/s)
Fd(X)=1/(1+Kd*X[14]); #DIP/TIP
FPNH4(X)=X[3]*X[4]/((KmN+X[3])*(KmN+X[4]))+X[3]*KmN/((X[3]+X[4])*(KmN+X[4])); #Preference between N compounds
fphs(X)=1/12*(32+48*(1-FPNH4(X))); #Photosyntesis stochiometry coefficient
VC(X)=(X[12]+X[13])/(KVC+X[12]+X[13]); #Vegetation covering
GL(X,t)=2.718*fp(t)*(1-VC(X))/((Ke0+aXp*X[7])*h(t))*(exp(-I0(t)/Is*exp(-(Ke0+aXp*X[7])*h(t)))-exp(-I0(t)/Is));
IGS(t)=max(sin(2*pi*t/tabs[end])*(sin(2*pi*t/tabs[end])>0),0.2); #Vegetation growth function
IGD(t)=max(-sin(2*pi*t/tabs[end])*(-sin(2*pi*t/tabs[end])>0),0.2); #Vegetation degradation function
DOsat(t)=0.475/(T(t)-239.65); #Dissolved oxygen saturation concentration (kg/m3)
Rw(t)=Ew*(tauw(t)/tauwcr(t)-1)*S(h(t))/V(h(t))*(tauw(t)>tauwcr(t)); #Resuspension by wind (kg/m3/s)
Xg(X)=X.*[Fd(X),0,1,1,1,1,0,1,1,0,0,0,0,0,0,1]; #Output concentrations to groundwater (kg/m3)
Xw(X)=X.*[1,0,1,1,1,1,1,1,1,1,1,0,0,1,1,1]; #Output concentrations (kg/m3)
psi1(X,t)=Kh*thetakh^(T(t)-293.15)*X[10]/(Kx*X[8]+X[10])*X[16]/(KhydDO+X[16])*X[8]; #Hydrolysis
psi2(X,t)=Kh*Oh*thetakh^(T(t)-293.15)*X[10]/(Kx*X[8]+X[10])*KhydDO/(KhydDO+X[16])*X[4]/(KNOh+X[4])*X[8]; #Anoxic hydrolysis
psi3(X,t)=muH*thetamuH^(T(t)-293.15)*X[5]/(Ks+X[5])*X[16]/(KOH+X[16])*X[3]/(KNHH+X[3])*Fd(X)*X[1]/(KPH+Fd(X)*X[1])*(1-(X[8]+X[9])/Kbiomax)*X[8]; #Aerobic growth of XH
psi4(X,t)=etaNO3*muH*thetamuH^(T(t)-293.15)*X[5]/(Ks+X[5])*KOH/(KOH+X[16])*X[4]/(KOH+X[4])*X[3]/(KNHH+X[3])*Fd(X)*X[1]/(KPH+Fd(X)*X[1])*(1-(X[8]+X[9])/Kbiomax)*X[8]; #Anoxic growth of XH
psi5(X,t)=bH*thetabH^(T(t)-293.15)*X[8]; #Lysis of XH
psi6(X,t)=muA*thetamuA^(T(t)-293.15)*X[3]/(KNHA+X[3])*X[16]/(KOA+X[16])*Fd(X)*X[1]/(KPA+Fd(X)*X[1])*(1-(X[8]+X[9])/Kbiomax)*X[9]; #Growth of XA
psi7(X,t)=bA*thetabA^(T(t)-293.15)*X[9]; #Lysis of XA
psi8(X,t)=1/h(t)*Kpl*thetaupl^(T(t)-293.15)*X[3]/(KNHP+X[3])*Fd(X)*X[1]/(KPP+Fd(X)*X[1])*(1-X[12]/KXml)*X[12]*IGS(t); #NH4 uptake by Xml
psi9(X,t)=1/h(t)*Kpl*thetaupl^(T(t)-293.15)*X[4]/(KNOP+X[4])*KNHP/(KNHP+X[3])*Fd(X)*X[1]/(KPP+Fd(X)*X[1] )*(1-X[12]/KXml)*X[12]*IGS(t); #NO3 uptake by Xml
psi10(X,t)=1/h(t)*Kdeg*thetadeg^(T(t)-293.15)*X[13]; #Xmd degradation
psi11(X,t)=S(h(t))/(0.1*V(h(t)))*DNH4*thetaNH4^(T(t)-293.15)*(NH4sed-X[3]); #Diffusion of NH4
psi12(X,t)=S(h(t))/(0.1*V(h(t)))*DNO3*thetaNO3^(T(t)-293.15)*(NO3sed-X[4]); #Diffusion of NO3
psi13(X,t)=vsx/h(t)*(1+Kvegsed*VC(X))*X[10]; #Sedimentation of Xs
psi14(X,t)=vsx/h(t)*(1+Kvegsed*VC(X))*X[11]; #Sedimentation of Xi
psi15(X,t)=Pmaxuptake*(Pmax-X[2])/(Pmax-Pmin)*(Fd(X)*X[1])/(kDIPup+Fd(X)*X[1]); #Phosphorus uptake by Xp
psi16(X,t)=Gmax*thetaG^(T(t)-293.15)*(X[2]-Pmin)/(Pmax-Pmin)*X[7]*GL(X,t); #Growth of Xp
psi17(X,t)=Kr*thetar^(T(t)-293.15)*X[7]; #Decay of Xp
psi18(X,t)=Kr*thetar^(T(t)-293.15)*X[2]; #Lysis of Pint
psi19(X,t)=Kresp*thetaresp^(T(t)-293.15)*X[16]/(KDOXp+X[16])*X[7]; #Respiration of Xp
psi20(X,t)=vsXp/h(t)*X[7]; #Sedimentation of Xp
psi21(X,t)=vsTIP/h(t)*(1-Fd(X))*(1+Kvegsed*VC(X))*X[1]; #Sedimentation of PIP
psi22(X,t)=DDIP*Kdif*thetadif^(T(t)-293.15)*(Psed-X[1]*Fd(X))*1/(0.1*h(t)); #Diffusion of DIP
psi23(X,t)=vsXTSS/h(t)*(1+Kvegsed*VC(X))*X[14]; #Sedimentation of vsXTSS
psi24(X,t)=Kavi*LinearAlgebra.dot(Ne,RF)/S(h(t))*(1-Kvegres*VC(X)); #Resuspension by avifauna
psi25(X,t)=KminOP*thetaminOP^(T(t)-293.15)*X[16]/(X[16]+kDO)*X[15]; #Mineralization of OP
psi26(X,t)=vsOP/h(t)*fPOP*(1+Kvegsed*VC(X))*X[15]; #Sedimentation of OP
psi27(X,t)=Rw(t)*(1-Kvegres*VC(X)); #Resuspension of TSS
psi28(X,t)=Rw(t)*(1-Kvegres*VC(X))*iPsed*Fpr; #Resuspension of OP
psi29(X,t)=Rw(t)*(1-Kvegres*VC(X))*iPsed*(1-Fpr); #Resuspension of PIP
psi30(X,t)=krea(t)*(DOsat(t)-X[16]); #Reaireation
psi31(X,t)=KphXp*thetaG^(T(t)-293.15)*X[7]; #Photosyntesis of Xp
psi32(X,t)=KphXml*thetaupl^(T(t)-293.15)*X[12]; #Photosyntesis of Xml
psi33(X,t)=Kmlmd*thetamlmd^(T(t)-293.15)*X[12]; #Degradation of Xml
#Processes
psi(X,t)=[psi1(X,t),psi2(X,t),psi3(X,t),psi4(X,t),psi5(X,t),psi6(X,t),psi7(X,t),psi8(X,t),psi9(X,t),psi10(X,t),psi11(X,t),psi12(X,t),psi13(X,t),psi14(X,t),psi15(X,t),psi16(X,t),psi17(X,t),psi18(X,t),psi19(X,t),psi20(X,t),psi21(X,t),psi22(X,t),psi23(X,t),psi24(X,t),psi25(X,t),psi26(X,t),psi27(X,t),psi28(X,t),psi29(X,t),psi30(X,t),psi31(X,t),psi32(X,t),psi33(X,t)];
#Stochiometry matrix
SM(X)=LinearAlgebra.Transpose([0 0 U1NH4 0 1-fhydSi fhydSi 0 0 0 -1 0 0 0 0 U1TIP -1;
U2TIP 0 U2NH4 0 1-fhydSi fhydSi 0 0 0 -1 0 0 0 0 0 0;
-iPbm+iPSs/YH 0 -iNbm+iNSs/YH 0 -1/YH 0 0 1 0 0 0 0 0 0 0 -1/YH+1;
-iPbm+iPSs/YH 0 -iNbm+iNSs/YH -(1-YH)/(2.86*YH) -1/YH 0 0 1 0 0 0 0 0 0 0 0;
U5TIP 0 U5NH4 0 fbmSs 0 0 -1 0 U5Xs fbmXi 0 0 0 0 0;
-iPbm 0 -iNbm-1/YA 1/YA 0 0 0 0 1 0 0 0 0 0 0 -4.57/YA+1;
U7TIP 0 U7NH4 0 fbmSs 0 0 0 -1 U7Xs fbmXi 0 0 0 0 0;
-iPXm 0 -iNXm 0 0 0 0 0 0 0 0 1 0 0 0 0;
-iPXm 0 0 -iNXm 0 0 0 0 0 0 0 1 0 0 0 0;
U10TIP 0 U10NH4 0 0 0 0 0 0 1-fpl fpl 0 -1 0 0 0;
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0;
-X[7] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 -iPXp/X[7] -iNXp*FPNH4(X) -iNXp*(1-FPNH4(X)) 0 0 1 0 0 0 0 0 0 iTSSXp 0 0;
U17TIP 0 U17NH4 0 fXpSs*iCODXp 0 -1 0 0 U17Xs fXpXi*iCODXp 0 0 -iTSSXp iPXp*Fop 0;
X[7] -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -iCODXp;
0 0 0 0 0 0 -1 0 0 0 0 0 0 -iTSSXp 0 0;
-1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0;
0 0 0 0 0 0 0 0 0 iCODsed 0 0 0 1 iPsed*Fpr 0;
1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 -fPDO;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
1-Fd(X) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 fphs(X);
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 fphs(X);
0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0]);
#Solver (BDF)
dX(X,p,t)=Qin(t)/V(h(t))*(Xin-Xw(X))+QET(h(t),t)/V(h(t))*Xw(X)-QP(t)/V(h(t))*Xw(X)-Qg(h(t))/V(h(t))*(Xg(X)-Xw(X))+SM(X)*psi(X,t);
@time X=DifferentialEquations.solve(DifferentialEquations.ODEProblem(dX,X0,(0,tabs[end])),DifferentialEquations.TRBDF2());
```

Any suggestions how could I optimize my code to be faster?

Thanks in advance for any advice!