ns: number of snapshots to store the solution
function numerical(nx,ny,nt,dx,dy,dt,re,wn,ns)
w1f = Array{Complex{Float64}}(undef,nx,ny)
w2f = Array{Complex{Float64}}(undef,nx,ny)
wnf = Array{Complex{Float64}}(undef,nx,ny)
j1f = Array{Complex{Float64}}(undef,nx,ny)
j2f = Array{Complex{Float64}}(undef,nx,ny)
jnf = Array{Complex{Float64}}(undef,nx,ny)
ut = Array{Float64}(undef, nx+1, ny+1)
data = Array{Complex{Float64}}(undef,nx,ny)
d1 = Array{Float64}(undef, nx, ny)
d2 = Array{Float64}(undef, nx, ny)
d3 = Array{Float64}(undef, nx, ny)
k2 = Array{Float64}(undef, nx, ny)
for i = 1:nx for j = 1:ny
data[i,j] = complex(wn[i+1,j+1],0.0)
end~end
wnf = fft(data)
k2 = wavespace(nx,ny,dx,dy) # wavespace over 2D domain
wnf[1,1] = 0.0
alpha1, alpha2, alpha3 = 8.0/15.0, 2.0/15.0, 1.0/3.0
gamma1, gamma2, gamma3 = 8.0/15.0, 5.0/12.0, 3.0/4.0
rho2, rho3 = -17.0/60.0, -5.0/12.0
for i = 1:nx for j = 1:ny
z = 0.5*dt*k2[i,j]/re
d1[i,j] = alpha1*z
d2[i,j] = alpha2*z
d3[i,j] = alpha3*z
end~end
for k = 1:nt
jnf = jacobian(nx,ny,dx,dy,wnf,k2)
for i = 1:nx for j = 1:ny
w1f[i,j] = ((1.0 - d1[i,j])/(1.0 + d1[i,j]))*wnf[i,j] +
(gamma1*dt*jnf[i,j])/(1.0 + d1[i,j])
end~end
w1f[1,1] = 0.0
j1f = jacobian(nx,ny,dx,dy,w1f,k2)
for i = 1:nx for j = 1:ny
w2f[i,j] = ((1.0 - d2[i,j])/(1.0 + d2[i,j]))*w1f[i,j] +
```

```
(rho2*dt*jnf[i,j] + gamma2*dt*j1f[i,j])/(1.0 + d2[i,j])
end~end
w2f[1,1] = 0.0
j2f = jacobian(nx,ny,dx,dy,w2f,k2)
for i = 1:nx for j = 1:ny
wnf[i,j] = ((1.0 - d3[i,j])/(1.0 + d3[i,j]))*w2f[i,j] +
(rho3*dt*j1f[i,j] + gamma3*dt*j2f[i,j])/(1.0 + d3[i,j])
end end
end
ut[1:nx,1:ny] = real(ifft(wnf)) # final solution
ut[nx+1,:] = ut[1,:]
ut[:,ny+1] = ut[:,1]
return ut
end