# Starting from MATLAB code, how do I use ode45 to solve Navier-Stokes in Julia?

**URL:** <https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187>\
**Category:** New to Julia\
**Tags:** question, diffeq, ode\
**Created:** [July 13, 2022, 4:56pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187 "2022-07-13T16:56:00Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [July 13, 2022, 4:56pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/1 "2022-07-13T16:56:00Z")

</div>

I’m a new Julia user but matlab expert(I am tired of paying them).

I want to call back certain values on my ode during solving. I have watched the tutorial from Chris [https://www.youtube.com/watch?v=KPEqYtEd-zY](https://www.youtube.com/watch?v=KPEqYtEd-zY) but not clear.

In matlab, I would do it as follows:

```julia
function [eqn,Qgap,Force,Fr]=navier(t,y)

%parameters
rho=810;
beta=5e8;
r1=0.01;
r2=0.004;
Ap=pi*(0.01^2-0.004^2);
L=0.06;
a=0.01;
f=10;
l=0.01;
Aeff=pi*2*r1*l;
delta=0.007;
N=20;
dx=delta/N;
rad=linspace(r1,r1+delta,N);

vel=y(1:N);
p1=y(N+1);
p2=y(N+2);

x=a*sin(2*pi*f*t);
u=a*2*pi*f*cos(2*pi*f*t);

eqn=zeros(N+2,1);
vel(1)=u;
vel(end)=0;
dudx=[(vel(2)-vel(1))/dx;(vel(3:end) -vel(1:end-2))/(2*dx);(vel(end)-vel(end-1))/dx];
dudx2=(vel(3:end)-2*vel(2:end-1)+vel(1:end-2))/(dx*dx);
dudx2 = [(dudx(1)-dudx(2))/dx;dudx2;(dudx(end-1)-dudx(end))/dx];
mu=[630 * (1+(0.05*(dudx)).^2).^(-0.23)];
pf=12*a*0.707*u*cos(2*pi*f*t)/(delta^2)*630*(1+(0.05*(0.707*u/delta)^2))^(-0.23);

Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta);
% Qgap=trapz(rad,2*pi*transpose(rad).*vel);

eqn(1)=(p1-p2+pf)/l + ((mu(1))/rho)*dudx2(1);
eqn(2:N-1)=(p1-p2+pf)/l + ((mu(2:N-1))./rho).*dudx2(2:N-1);
eqn(N)=(p1-p2+pf)/l + ((mu(N))/rho)*dudx2(N);
eqn(N+1)=(beta/(Ap*(L-l-x)))*(-Qgap+Ap*u);
eqn(N+2)=(beta/(Ap*(L-l+x)))*(Qgap-Ap*u);

Force=(p1-p2+pf)*Ap;
Fr=mu(1)*Aeff*dudx(1);

    
end

```

In matlab, we just have to put the variables as outputs after the eqn, then write a callback function (a bit complicated version of my own)

```julia
[~,Qgap,Force,Fr] = cellfun(@(tsol,ysol) navier(tsol,ysol.'), num2cell(tsol), num2cell(ysol,2),'uni',0);
Qgap=cell2mat(Qgap);
Force=cell2mat(Force);
Fr=cell2mat(Fr);

```

where as it gives me back a cell of the values with the time points, and then convert to a matlab simple array.

Can you help me achieve this in Julia?

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [July 13, 2022, 5:47pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/2 "2022-07-13T17:47:24Z")

</div>

You’re just trying to save a few extra variables of interest, right? You should be able to handle that with a [`SavingCallback`](https://diffeq.sciml.ai/stable/features/callback_library/#saving_callback). Do you have a working Julia version of the code without callbacks?

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [July 13, 2022, 7:32pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/3 "2022-07-13T19:32:43Z")

</div>

No I don’t have it. I was looking for every option of conversion from Matlab to Julia so I can do it smoothly.  
If you have got free time and alot of coding interest, please fell free.

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [July 13, 2022, 8:29pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/4 "2022-07-13T20:29:20Z")

</div>

It would be easier to help you if you could attempt some of the conversion yourself. For example, I think you could have ported `navier` for us and maybe asked how to do the numerical quadrature.

Nonetheless, I’m feeling generous. I’m attempted to port `navier`. I did not think too much about it, so there is likely an error:

```julia
using NumericalIntegration

function navier(t,y)

#parameters
rho=810
beta=5e8
r1=0.01
r2=0.004
Ap=pi*(0.01^2-0.004^2)
L=0.06
a=0.01
f=10
l=0.01
Aeff=pi*2*r1*l
delta=0.007
N=20
dx=delta/N
rad=range(r1,r1+delta,length=N)

vel=y[1:N]
p1=y[N+1]
p2=y[N+2]

x=a*sin(2*pi*f*t)
u=a*2*pi*f*cos(2*pi*f*t)

eqn=zeros(N+2,1)
vel[1]=u
vel[end]=0
dudx=[(vel[2]-vel[1])/dx;(vel[3:end] -vel[1:end-2])/(2*dx);(vel[end]-vel[end-1])/dx]
dudx2=(vel[3:end]-2*vel[2:end-1].+vel[1:end-2])/(dx*dx)
dudx2 = [(dudx[1]-dudx[2])/dx;dudx2;(dudx[end-1]-dudx[end])/dx]
mu=630 * (1 .+ (0.05*(dudx)).^2).^(-0.23)
pf=12*a*0.707*u*cos(2*pi*f*t)/(delta^2)*630*(1+(0.05*(0.707*u/delta)^2))^(-0.23)

#Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta)
# Qgap=trapz(rad,2*pi*transpose(rad).*vel)
Qgap = integrate(rad, 2*pi*rad.*vel)

eqn[1]=(p1-p2+pf)/l .+ ((mu[1])/rho)*dudx2[1]
eqn[2:N-1]=(p1-p2+pf)/l .+ ((mu[2:N-1])./rho).*dudx2[2:N-1]
eqn[N]=(p1-p2+pf)/l .+ ((mu[N])/rho)*dudx2[N]
eqn[N+1]=(beta/(Ap*(L-l-x)))*(-Qgap+Ap*u)
eqn[N+2]=(beta/(Ap*(L-l+x)))*(Qgap-Ap*u)

Force=(p1-p2+pf)*Ap
Fr=mu[1]*Aeff*dudx[1]

return (
    eqn=eqn,
    Qgap=Qgap,
    Force=Force,
    Fr=Fr
)

end

```

It would probably be best to use a for loop to aggregate the results, but I’ll show you how to do with `map`.

```julia
julia> results = map((t,y)->navier(t,y)[(:Qgap, :Force, :Fr)], 1:10, eachcol(rand(22,10)))
10-element Vector{NamedTuple{(:Qgap, :Force, :Fr), Tuple{Float64, Float64, Float64}}}:
 (Qgap = 0.00030102609426749455, Force = 53.32666898640483, Fr = 24.619494874431233)
 (Qgap = 0.00029967156523574053, Force = 53.32663257434115, Fr = -49.498430150318434)
 (Qgap = 0.0002801507212840063, Force = 53.32680302845422, Fr = 61.985357870666995)
 (Qgap = 0.00033119472944960455, Force = 53.327010098387575, Fr = 37.68745577297903)
 (Qgap = 0.00029727868343977574, Force = 53.326764122592515, Fr = -50.54111347534392)
 (Qgap = 0.0003496617278965436, Force = 53.3268011948002, Fr = 8.655337852075665)
 (Qgap = 0.0003150886996319404, Force = 53.32660908124239, Fr = -81.82836381327539)
 (Qgap = 0.00028876439037434895, Force = 53.32686035546042, Fr = 61.154294238471735)
 (Qgap = 0.0002664868528073511, Force = 53.326850589478646, Fr = -68.03440757714097)
 (Qgap = 0.00024699380480928317, Force = 53.326924487620666, Fr = 5.137304480117839)

julia> Qgap = getproperty.(results, :Qgap)
10-element Vector{Float64}:
 0.00030102609426749455
 0.00029967156523574053
 0.0002801507212840063
 0.00033119472944960455
 0.00029727868343977574
 0.0003496617278965436
 0.0003150886996319404
 0.00028876439037434895
 0.0002664868528073511
 0.00024699380480928317

julia> Force = getproperty.(results, :Force)
10-element Vector{Float64}:
 53.32666898640483
 53.32663257434115
 53.32680302845422
 53.327010098387575
 53.326764122592515
 53.3268011948002
 53.32660908124239
 53.32686035546042
 53.326850589478646
 53.326924487620666

julia> Fr = getproperty.(results, :Fr)
10-element Vector{Float64}:
  24.619494874431233
 -49.498430150318434
  61.985357870666995
  37.68745577297903
 -50.54111347534392
   8.655337852075665
 -81.82836381327539
  61.154294238471735
 -68.03440757714097
   5.137304480117839

```

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [July 13, 2022, 8:50pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/5 "2022-07-13T20:50:36Z")

</div>

Thank you very much for your generosity. I was looking forward to write it tomorrow morning, fresh start with a cup of coffee (Numerical Modelling is Numerically tiresome).

That’s all. No ODEProblem(navier,u0,tspan…), no choosing methods to solve?

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [July 13, 2022, 9:41pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/6 "2022-07-13T21:41:29Z")

</div>

I only parrotted back the code you gave us. If you want more, we’ll need more _executable_ MATLAB code.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [July 13, 2022, 9:51pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/7 "2022-07-13T21:51:01Z")

</div>

Here is the full code. I have altered quite few parameteres (prototype of the lab).

```julia
f=10;
a=0.01;
time=linspace(0.25/f,2/f,100);

% 
xglob=a*sin(2*pi*f*time);
uglob=a*2*pi*f*cos(2*pi*f*time);

N=100;
y0(1:N,1)=0;
y0(N+1)=8e5;
y0(N+2)=8e5;
tic
[tsol,ysol]=ode45(@navier,time,y0);
toc

[~,Qgap,Force,Fr] = cellfun(@(tsol,ysol) navier(tsol,ysol.'), num2cell(tsol), num2cell(ysol,2),'uni',0);
Qgap=cell2mat(Qgap);

figure(1)
hold on
yyaxis left
title('Qgap and Velocity')
xlabel('time')
ylabel('V(m/s)')
plot(tsol,ysol(:,1),'b--',tsol,ysol(:,N/2),'r--',tsol,ysol(:,N),'m--',tsol,uglob);            
yyaxis right
ylabel('m^3/s')
plot(tsol,Qgap,'ko-');

hold off

Force=cell2mat(Force);
Fr=cell2mat(Fr);

figure(2)
hold on
yyaxis left
title('Chamber Pressures')
xlabel('time[s]')
ylabel('Pressure[Pa]')
plot(tsol,ysol(:,N+1),'k',tsol,ysol(:,N+2));            
yyaxis right
ylabel('dist[m]')
legend('ComCha','RebCha')
set(gcf,'Position', [2 2 1000 1000])
plot(tsol,xglob);

hold off

function [eqn,Qgap,Force,Fr]=navier(t,y)

%parameters
rho=810;
beta=5e8;
r1=0.01;
r2=0.004;
Ap=pi*(0.01^2-0.004^2);
L=0.06;
a=0.01;
f=10;
l=0.01;
Aeff=pi*2*r1*l;
delta=0.007;
N=100;
dx=delta/N;
rad=linspace(r1,r1+delta,N);

vel=y(1:N);
p1=y(N+1);
p2=y(N+2);

x=a*sin(2*pi*f*t);
u=a*2*pi*f*cos(2*pi*f*t);

eqn=zeros(N+2,1);

dudx=[(u-vel(1))/dx;(vel(3:end) -vel(1:end-2))/(2*dx);(vel(end)-0)/dx];
dudx2=(vel(3:end)-2*vel(2:end-1)+vel(1:end-2))/(dx*dx);
dudx2 = [(dudx(1)-dudx(2))/dx;dudx2;(dudx(end-1)-dudx(end))/dx];
mu=[630 * (1+((0.05*(dudx)).^2)).^(-0.23)];
pf=12*a*0.707*u*cos(2*pi*f*t)/(delta^2)*630*(1+(0.05*(0.707*u/delta)^2))^(-0.23);

eqn(1)=(p1-p2+pf)/l + ((mu(1))/rho)*dudx2(1);
eqn(2:N-1)=(p1-p2+pf)/l + ((mu(2:N-1))./rho).*dudx2(2:N-1);
eqn(N)=(p1-p2+pf)/l + ((mu(N))/rho)*dudx2(N);

Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta);
% Qgap=trapz(rad,2*pi*transpose(rad).*vel);

eqn(N+1)=(beta/(Ap*(L-x)))*(-Qgap+Ap*u);
eqn(N+2)=(beta/(Ap*(L+x)))*(Qgap-Ap*u);

Force=(p1-p2+pf)*Ap;
Fr=mu(1)*Aeff*dudx(1);

%eqn = transpose([temp(1,:) temp(2,:) temp(3,:)]);
    

end

```

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [July 14, 2022, 8:11am UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/8 "2022-07-14T08:11:05Z")

</div>

I suggest posting the example as a new post asking specifically about ODE solvers. Perhaps you could even combine the video from Chris with the Julia function I ported for you to make an attempt and then ask a more specific question.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [July 15, 2022, 7:01pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/9 "2022-07-15T19:01:50Z")

</div>

That makes it more complicated don’t you think. I wish we worked on this, that it makes it easier for the other users new and experts out there.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [July 15, 2022, 7:14pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/10 "2022-07-15T19:14:25Z")

</div>

If you want to learn how to solve ODEs with callbacks in DifferentialEquations.jl I suggest carefully going through the docs at [Event Handling and Callback Functions · SciML](https://docs.sciml.ai/dev/modules/DiffEqDocs/features/callback_functions/) and at [Callback Library · SciML](https://docs.sciml.ai/dev/modules/DiffEqDocs/features/callback_library/). These have examples showing how to solve differential equations using various kinds of callbacks.

Once you’ve read them, if you want to put together a _simple_ illustrative Julia example and post it here as a new issue I’m sure you’ll be able to get help with any problems that you may encounter solving ODEs in Julia. Getting a simple example running with the features you need (i.e. callbacks and such) would allow you to understand the DifferentialEquations.jl workflow, and presumably make it easier for you to port your more complicated codes over.

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [July 15, 2022, 11:11pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/11 "2022-07-15T23:11:46Z")

</div>

> [@simonluzara](#):
>
> That makes it more complicated don’t you think. I wish we worked on this, that it makes it easier for the other users new and experts out there.

I changed the title. Maybe @ChrisRackauckas or someone else is feeling generous enough to port the code further for you.

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [July 15, 2022, 11:19pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/12 "2022-07-15T23:19:43Z")

</div>

Per [https://diffeq.sciml.ai/stable/solvers/ode\_solve/](https://diffeq.sciml.ai/stable/solvers/ode_solve/)

> For users familiar with MATLAB/Python/R, good translations of the standard library methods are as follows:
> 
> - `ode45`/`dopri5` –\> `DP5()`, though in most cases `Tsit5()` is more efficient

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 16, 2022, 5:07am UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/13 "2022-07-16T05:07:58Z")

</div>

Just make it `navier(dy,y,p,t)` from @mkitti 's code and you’re done. Did you try that?

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [September 7, 2022, 12:44pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/15 "2022-09-07T12:44:00Z")

</div>

Hi Chris. Been tough with company meetings and holidays. Came back fresh and attacked the problem. But still I have got a bug.

I have followed your tutorial in setting up constant parameters and calculated ones inside the function from [https://tutorials.sciml.ai/html/introduction/03-optimizing\_diffeq\_code.html](https://tutorials.sciml.ai/html/introduction/03-optimizing_diffeq_code.html)

```julia
p = (1.0,1.0,1.0,10.0,0.001,100.0,Ayu,uAx,Du,Ayv,vAx,Dv) # a,α,ubar,β,D1,D2
function gm4!(dr,r,p,t)
  a,α,ubar,β,D1,D2,Ayu,uAx,Du,Ayv,vAx,Dv = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  mul!(Ayu,Ay,u)
  mul!(uAx,u,Ax)
  mul!(Ayv,Ay,v)
  mul!(vAx,v,Ax)
  @. Du = D1*(Ayu + uAx)
  @. Dv = D2*(Ayv + vAx)
  @. du = Du + a*u*u./v + ubar - α*u
  @. dv = Dv + a*u*u - β*v
end
prob = ODEProblem(gm4!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

```

So I have done the same with my code, but I get an error;

```julia
#parameters
p=(810,5e8,0.01,0.004,pi*(0.01^2-0.004^2),0.06,0.01,10,0.01,pi*2*0.01*0.01,0.007,100,0.0007/100,
    range(0.01,0.017,length=100),x,u,μ,Pf,Qg)

function navier!(dy,y,p,t)
    #parameters
    ρ,β,r₁,r₂,Aₚ,L₀,amp,f,Lₚ,Aₑ,δ,N,Δx,rad,x,u,μ,Pf,Qg=p
    vel=y[1:N]
    p1=y[N+1]
    p2=y[N+2]

    x=amp*sin(2*pi*f*t)
    u=amp*2*pi*f*cos(2*pi*f*t)

    dudx=[(u-vel(1))/Δx;(vel(3:end) -vel(1:end-2))/(2*Δx);(vel(end)-0)/Δx]
    dudx2=(vel[3:end]-2*vel[2:end-1].+vel[1:end-2])/(Δx*Δx)
    dudx2 = [(dudx[1]-dudx[2])/Δx;dudx2;(dudx[end-1]-dudx[end])/Δx]
    μ=630 * (1 .+ (0.05*(dudx)).^2).^(-0.23)
    Pf=12*amp*0.707*u*cos(2*pi*f*t)/(δ^2)*630*(1+(0.05*(0.707*u/δ)^2))^(-0.23)

    #Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta)
    # Qgap=trapz(rad,2*pi*transpose(rad).*vel)
    Qgap = simps(2*pi*rad.*vel,r₁,r₁+δ)

    dy[1:end-2]=(p1-p2+pf)/Lₚ. + ((μ[1:N])./rho).*dudx2[1:N]
    dy[end-1]=(β/(Aₚ*(L₀-Lₚ/2-x)))*(-Qgap+Aₚ*u)
    dy[end]=(β/(Aₚ*(L₀-Lₚ/2+x)))*(Qgap-Aₚ*u)

    Force=(p1-p2+Pf)*Aₚ
    Fr=μ[1]*Aₑ*dudx[1]

    return (Qgap=Qgap,Force=Force,Fr=Fr)

end

```

> UndefVarError: x not defined
> 
> Stacktrace:  
> [1] top-level scope  
> @ In[9]:2  
> [2] eval  
> @ .\boot.jl:373 [inlined]  
> [3] include\_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)  
> @ Base .\loading.jl:1196

Can you or anyone help me please?

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [September 7, 2022, 1:13pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/16 "2022-09-07T13:13:19Z")

</div>

I haven’t read this thread and am not sure whether the code snippet you posted is supposed to be a self contained MWE, but if it is the error is obvious:

```julia
#parameters
p=(810,5e8,0.01,0.004,pi*(0.01^2-0.004^2),0.06,0.01,10,0.01,pi*2*0.01*0.01,0.007,100,0.0007/100,
    range(0.01,0.017,length=100),x,u,μ,Pf,Qg)

```

here you are referencing a variable `x` which is apparently undefined. In the notebook you linked, `Ayu`, `uAx` etc are all defined in previous cells, so you would similarly have to define `x`, `u`, `μ`, `Pf`, and `Qg` before using them.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [September 7, 2022, 1:37pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/17 "2022-09-07T13:37:28Z")

</div>

Yes I thought of that too. Since x is dependent in time, it will be an iterative parameter. I have removed it from the list of p with u,Pf,Qg.

Now I get another error which I think comes from me not knowing to create non-allocation arrays inside the function. Here is the modified code.

```julia
#parameters
p=(810,5e8,0.01,0.004,pi*(0.01^2-0.004^2),0.06,0.01,10,0.01,pi*2*0.01*0.01,0.007,100,0.007/100,0.01:0.007/100:0.017)

function navier!(dy,y,p,t)
    #parameters
    ρ,β,r₁,r₂,Aₚ,L₀,amp,f,Lₚ,Aₑ,δ,N,Δx,rad=p
    vel=y[1:N]
    p1=y[N+1]
    p2=y[N+2]

    x=amp*sin(2*pi*f*t)
    u=amp*2*pi*f*cos(2*pi*f*t)

    dudx=[(u-vel(1))/Δx;(vel(3:end)-vel(1:end-2))/(2*Δx);(vel(end)-0)/Δx]
    dudx2=(vel[3:end]-2*vel[2:end-1].+vel[1:end-2])/(Δx*Δx)
    dudx2 = [(dudx[1]-dudx[2])/Δx;dudx2;(dudx[end-1]-dudx[end])/Δx]
    μ=630 * (1 .+ (0.05*(dudx)).^2).^(-0.23)
    Pf=12*amp*0.707*u*cos(2*pi*f*t)/(δ^2)*630*(1+(0.05*(0.707*u/δ)^2))^(-0.23)

    #Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta)
    # Qgap=trapz(rad,2*pi*transpose(rad).*vel)
    Qgap = simps(2*pi*rad.*vel,r₁,r₁+δ)

    dy[1:end-2]=(p1-p2+pf)/Lₚ. + ((μ[1:N])./rho).*dudx2[1:N]
    dy[end-1]=(β/(Aₚ*(L₀-Lₚ/2-x)))*(-Qgap+Aₚ*u)
    dy[end]=(β/(Aₚ*(L₀-Lₚ/2+x)))*(Qgap-Aₚ*u)

    Force=(p1-p2+Pf)*Aₚ
    Fr=μ[1]*Aₑ*dudx[1]

    return (Qgap=Qgap,Force=Force,Fr=Fr)

end

```

I am getting this error;

> syntax: missing last argument in “3:” range expression
> 
> Stacktrace:  
> [1] top-level scope  
> @ In[3]:14  
> [2] eval  
> @ .\boot.jl:373 [inlined]  
> [3] include\_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)  
> @ Base .\loading.jl:1196

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [September 7, 2022, 1:42pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/18 "2022-09-07T13:42:53Z")

</div>

> [@simonluzara](#):
>
> `vel(3:end)`

Julia uses square brackets for array indexes [and] go through and change all instances of using ( and )

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [September 7, 2022, 1:53pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/19 "2022-09-07T13:53:06Z")

</div>

I feel dumb as a primary school kid not seeing this 😀.

I think now the generic function has been created with this code.  
Two things;

1. I want later on to get to the origin of this post (calling the variables Qgap,Force and Fr after the solution).
2. Any suggestion to optimise the code, as I want to increase the FDM spacing to N=1000 for stability reasons and am afraid it will take ages.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [September 7, 2022, 3:04pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/20 "2022-09-07T15:04:27Z")

</div>

the main thing is to stop rolling your own integration. [Solving Integro Differential Equations · NeuralPDE.jl](https://neuralpde.sciml.ai/stable/tutorials/integro_diff/) describes how to deal with integrodifferential equations more efficiently.

---

<div class="post-metadata">

**Author:** ![simonluzara](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simonluzara/32/205260_2.png) [@simonluzara](https://discourse.julialang.org/u/simonluzara)\
**Post date:** [September 7, 2022, 3:17pm UTC](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187/21 "2022-09-07T15:17:46Z")

</div>

Please be generous enough to help me through it. Here is the finest version.

```julia
#parameters
N=101
p=(810,5e8,0.01,0.004,pi*(0.01^2-0.004^2),0.06,0.01,10,0.01,pi*2*0.01*0.01,0.007,N,0.007/N,range(0.01,0.017,length=N))

function navier!(dy,y,p,t)
    #parameters
    ρ,β,r₁,r₂,Aₚ,L₀,amp,f,Lₚ,Aₑ,δ,N,Δx,rad=p
    vel=y[1:N]
    p1=y[N+1]
    p2=y[N+2]

    x=amp*sin(2*pi*f*t)
    u=amp*2*pi*f*cos(2*pi*f*t)

    dudx=[(u-vel[1])/Δx;(vel[3:end]-vel[1:end-2])/(2*Δx);(vel[end]-0)/Δx]
    dudx2=(vel[3:end]-2*vel[2:end-1].+vel[1:end-2])/(Δx*Δx)
    dudx2 = [(dudx[1]-dudx[2])/Δx;dudx2;(dudx[end-1]-dudx[end])/Δx]
    μ=630 * (1 .+ (0.05*(dudx)).^2).^(-0.23)
    Pf=12*amp*0.707*u*cos(2*pi*f*t)/(δ^2)*630*(1+(0.05*(0.707*u/δ)^2))^(-0.23)

    #Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta)
    # Qgap=trapz(rad,2*pi*transpose(rad).*vel)
    Qgap = simps(2*pi*rad.*vel,r₁,r₁+δ)

    dy[1:N]=(p1-p2+Pf)/Lₚ.+((μ[1:N])./ρ).*dudx2[1:N]
    dy[N+1]=(β/(Aₚ*(L₀-Lₚ/2-x)))*(-Qgap+Aₚ*u)
    dy[N+2]=(β/(Aₚ*(L₀-Lₚ/2+x)))*(Qgap-Aₚ*u)

    Force=(p1-p2+Pf)*Aₚ
    Fr=μ[1]*Aₑ*dudx[1]

    return (Qgap=Qgap,Force=Force,Fr=Fr)

end

```

with the solving method chosen as follows;

```julia
y0=zeros(Float64,N+2)
y0[1:N].=0
y0[N+1]=8e5
y0[N+2]=8e5
time=(0.0,0.5)
ode_prob=ODEProblem(navier!,y0,time,p);
@btime soln=solve(ode_prob,Rodas5(),abstol=1e-6,reltol=1e-6,saveat=0.3:0.001:0.4);

```

[Next page](https://discourse.julialang.org/t/starting-from-matlab-code-how-do-i-use-ode45-to-solve-navier-stokes-in-julia/84187.md?page=2)