# Flow driven by an oscillating solid cylinder and bounded by a stationary outer hollow cylinder
\begin{equation}\label{p1introhpa} \Psi _{\eta \eta \tau }^{\mathrm{Sym}}  = \frac{1}{{R_s }}\Psi^{\mathrm{Sym}}_{\eta \eta \eta \eta }  + \Psi^{\mathrm{Sym}}_{\eta \eta \eta } \Psi _\theta^{\mathrm{Sym}}   - \Psi^{\mathrm{Sym}}_{\eta \eta \theta } \Psi^{\mathrm{Sym}}_\eta, \\
\quad \Psi^{\mathrm{Sym}}(\eta= 0) = 0,\quad \Psi^{\mathrm{Sym}}_{\eta\eta}(\eta= 0) = 0, \qquad \Psi^{\mathrm{Sym}}(\eta= 1) = 0, \quad \Psi^{\mathrm{Sym}}_\eta(\eta= 1) = - \mu sin\theta, \qquad
0 \le \eta \le1, \quad 0 \le \theta \le \pi 
\end{equation}

## Calculating $\Psi(\theta,\eta)$ the steady state solution $\theta$-marching
The governing equation is solved using: 

* Half domain $\eta$-domain $0 \le \eta \le 1$
* $\theta$ marching with some of the nonlinear terms being explicit 


## Summary of results
The numerical soltion is only smooth for  $0\leq\theta\leq\pi$ for very small reynolds numbers $-1< R_s< 2.0$, at higher Reynolds numbers the numerical errors accumulated as $\theta$ increases from zero and seem to be related to the choice of the meshsize $\Delta\theta$. When $\Delta\theta$ is *increased* the errors seem to vanish.

The solution generated by $\theta$-marching still matches the time marched solution near $\theta = 0$ upto  around $R_s = 45$ with very strong agreement with the 1st Order winding timemarched solution and less so with the second order winding scheme. 
        
## Small $R_s$ series solution



```julia
using DifferentialEquations
using Plots
using BenchmarkTools
using DelimitedFiles
using LaTeXStrings
using SparseArrays
using LinearAlgebra
using StaticArrays
using Crayons.Box
using Printf
```


```julia
function ψ_ηηηη(n ::Int64, m:: Int64, ψ::Matrix{Float64}, meshsize::SVector{2, Int64}, delta_sizes::SVector{2, Float64}, paramtrs::Tuple{Float64, String})
     μ = paramtrs[1]
     order = paramtrs[2]
     N = meshsize[1]
    (Δθ, Δη) = delta_sizes
    
    if  n > 2 && n < N 
        
        return( ψ[n-2,m]  -4*ψ[n-1,m] +6*ψ[n,m] -4*ψ[n+1,m] + ψ[n+2,m])
        
    elseif n == 2 #Second order only
        
        return( (-ψ[n,m]) -4*ψ[n-1,m] +6*ψ[n,m] -4*ψ[n+1,m] + ψ[n+2,m])
    
    elseif n == N #Second order only
        
        return( ψ[n-2,m]  -4*ψ[n-1,m] +6*ψ[n,m] -4*ψ[n+1,m] + (ψ[n,m] - 2.0*Δη*μ*sin(Δθ*(m-1))))
    else
        println("Error in ψ_ηηη n =$n")
        return 0.0
    end
end



function ψ_ηηη(n ::Int64, m:: Int64, ψ::Matrix{Float64}, meshsize::SVector{2, Int64}, delta_sizes::SVector{2, Float64}, paramtrs::Tuple{Float64, String})
     μ = paramtrs[1]
     order = paramtrs[2]
     N = meshsize[1]
    (Δθ, Δη) = delta_sizes
    
    if  n > 2 && n < N 
        if n > 3 && n < N - 1 && order == "3rdOrder"
          return( ψ[n-3,m]/8 - ψ[n-2,m] + 13*ψ[n-1,m]/8 - 13*ψ[n+1,m]/8 + ψ[n+2,m] - ψ[n+3,m]/8)  
        else
          return(-0.5*ψ[n-2,m] + ψ[n-1,m] - ψ[n+1,m] + 0.5*ψ[n+2,m])
        end
        
    elseif n == 2
        
        return(-0.5*(-ψ[n,m]) + ψ[n-1,m] - ψ[n+1,m] + 0.5*ψ[n+2,m])
    
    elseif n == N
        
        return(-0.5*ψ[n-2,m] + ψ[n-1,m] - ψ[n+1,m] + 0.5*(ψ[n,m] - 2.0*Δη*μ*sin(Δθ*(m-1))))
    else
        println("Error in ψ_ηηη n =$n")
        return 0.0
    end
end



function ψ_ηη(n ::Int64, m:: Int64, ψ::Matrix{Float64}, meshsize::SVector{2, Int64})
    N = meshsize[1]
    if  n == 2 || n == N  
            
        return(  ψ[n-1,m] -2*ψ[n,m]+ ψ[n+1,m] )
    elseif n > 2  && n < N 
        return ( -ψ[n-2,m]/12 +4*ψ[n-1,m]/3 -2.5*ψ[n,m]+ 4*ψ[n+1,m]/3 -ψ[n+2,m]/12)
    else
        println("Error in ψ_ηη n =$n")
        return 0.0
    end
end


function ψ_η(n ::Int64, m:: Int64, ψ::Matrix{Float64},  meshsize::SVector{2, Int64})
    N = meshsize[1]
    if  n == 2 || n == N
        return(  -0.5*ψ[n-1,m] + 0.5*ψ[n+1,m] )
        
    elseif  n>=3 && n<= N-1
        return(  ψ[n-2,m]/12 -2*ψ[n-1,m]/3 + 2*ψ[n+1,m]/3 - ψ[n+2,m]/12 )
        
    else
        println("Error in ψ_η n =$n")
        return 0.0
    end
end

"""The symmetric (odd in η) asymptotic solution for the streamfunction for small \$R_s\$ with 5 terms"""    
function Ψ_smallr(θ::Float64, η::Float64 ,r::Float64, number_of_terms::Int) 
    terms = [
        -η*(η^2-1)*sin(θ)/2.0,
        
        r*η*(η^2-1)^2.0*(2.0+η^2)*sin(2*θ)/560.0, 
        
        r^2.0*η*(η^2-1)^2*((-591 -2294*η^2 + 161*η^4 + 1428*η^6)*sin(θ) + 
            +3.0*(423 + 166*η^2 - 553*η^4 - 84*η^6)sin(3*θ))/62092800.0,
        
        r^3.0*η*(η - 1)^2*(η + 1)^2*(2184*η^10*cos(θ)^2 + 1540*η^10 + 24213*η^8*cos(θ)^2 - 44205*η^8 + 
            + 18942*η^6*cos(θ)^2 + 33446*η^6 - 101444*η^4*cos(θ)^2 + 66052*η^4 - 23086*η^2*cos(θ)^2 - 
            6746*η^2 - 14889*cos(θ)^2 - 9383)*sin(θ)*cos(θ)/28252224000,
        
        r^4.0*η*(η - 1)^2*(η + 1)^2*(968647680*η^14*sin(θ)^5 + 14774464320*η^14*sin(θ) + 3510571680*η^14*sin(3*θ) -
            393513120*η^14*sin(5*θ) - 39671478000*η^12*sin(θ)^5 - 75537556710*η^12*sin(θ) - 
            28276357725*η^12*sin(3*θ) - 3752563815*η^12*sin(5*θ) - 181021246560*η^10*sin(θ)^5 + 
            357286638324*η^10*sin(θ) + 75142631718*η^10*sin(3*θ) - 1180776366*η^10*sin(5*θ) + 
            677920708080*η^8*sin(θ)^5 - 645834258762*η^8*sin(θ) - 123813530379*η^8*sin(3*θ) - 
            1711397457*η^8*sin(5*θ) + 193167076032*η^6*sin(θ)^5 - 82870720296*η^6*sin(θ) + 
            215941919508*η^6*sin(3*θ) + 59820477948*η^6*sin(5*θ) + 153978403824*η^4*sin(θ)^5 - 
            123221433370*η^4*sin(θ) + 181687022685*η^4*sin(3*θ) - 60662329497*η^4*sin(5*θ) - 
            500785425504*η^2*sin(θ)^5 + 367229681332*η^2*sin(θ) - 134582502330*η^2*sin(3*θ) - 
            33288605262*η^2*sin(5*θ) - 1507832445168*sin(θ)^5 + 969877856394*sin(θ) - 482955016725*sin(3*θ) - 
            9092499471*sin(5*θ))/16189310893916160000
        ]
    if number_of_terms < 6
        return sum(terms[1:number_of_terms])
    else
        return term[1]    
    end
end

"""First η-derivative of the symmetric streamfunction calculated from the small \$R_s\$ asymptotic series (3 terms)"""
function Ψ_η_smallr(θ::Float64, η::Float64 ,r::Float64) 
    return(-(3.0*η^2-1)*sin(θ)/2.0 + 
        +r*(7*η^6-9*η^2+2)*sin(2*θ)/560.0 + 
        + r^2.0*(3*(-197 - 1112*η^2 + 6930*η^4 - 2772*η^6 - 8085*η^8 + 5236*η^10)*sin(θ) + (1269 - 6120*η^2 - 6930*η^4 + 24948*η^6 - 10395*η^8 - 2772*η^10)*sin(3*θ))/62092800.0)
end

"""Second η-derivative of the symmetric streamfunction calculated from the small \$R_s\$ asymptotic series (3 terms)"""        
function  Ψ_ηη_smallr(θ::Float64, η::Float64 ,r::Float64) 
    return(-3.0*η*sin(θ) + r*(42*η^5-18*η)*sin(2*θ)/560.0 + 
        r^2.0*(3*( - 2224*η + 27720*η^3 - 16632*η^5 - 64680*η^7 + 52360*η^9)*sin(θ) + 
            (- 12240*η - 27720*η^3 + 149688*η^5 - 83160*η^7 - 27720*η^9)*sin(3*θ))/62092800.0)
end        

"""Third η-derivative of the symmetric streamfunction calculated from the small \$R_s\$ asymptotic series (3 terms)"""
function  Ψ_ηηη_smallr(θ::Float64, η::Float64 ,r::Float64) 
    return(-3.0*sin(θ) + r*(210*η^4-18)*sin(2*θ)/560.0 + 
        r^2.0*(3*( - 2224 + 83160*η^2 - 83160*η^4 - 452760*η^6 + 471240*η^8)*sin(θ) + (- 12240 - 83160*η^2 + 748440*η^4 - 582120*η^6 - 249480*η^8)*sin(3*θ))/62092800.0)
end 

"""Fourth η-derivative of the symmetric streamfunction calculated from the small \$R_s\$ asymptotic series (3 terms)"""
function  Ψ_ηηηη_smallr(θ::Float64, η::Float64 ,r::Float64) 
    return(3*r*η*(3*r*η^6*sin(θ)^2/70 + r*η^6/35 + 3*r*η^4*sin(θ)^2/40 - r*η^4/10 - 9*r*η^2*sin(θ)^2/140 + 3*r*η^2/70 + r*sin(θ)^2/280 + η^2*cos(θ))*sin(θ))
end 

"""1st θ-derivative of the symmetric streamfunction calculated from the small \$R_s\$ asymptotic series (3 terms)"""
function  Ψ_θ_smallr(θ::Float64, η::Float64 ,r::Float64)  
return (-η*(η - 1)*(η + 1)*(756*r^2*η^8*cos(θ)^3 - 924*r^2*η^8*cos(θ) + 4221*r^2*η^6*cos(θ)^3 - 2849*r^2*η^6*cos(θ) - 6471*r^2*η^4*cos(θ)^3 + 5467*r^2*η^4*cos(θ) - 2313*r^2*η^2*cos(θ)^3 + 1309*r^2*η^2*cos(θ) + 3807*r^2*cos(θ)^3 - 3003*r^2*cos(θ) - 55440*r*η^4*cos(2*θ) - 55440*r*η^2*cos(2*θ) + 110880*r*cos(2*θ) + 7761600*cos(θ))/15523200)
end    

"""ηηθ-derivative of the symmetric streamfunction calculated from the small \$R_s\$ asymptotic series (3 terms)"""
function  Ψ_ηηθ_smallr(θ::Float64, η::Float64 ,r::Float64)  
return (η*(3465*r^2*η^8*sin(θ)^2*cos(θ) + 770*r^2*η^8*cos(θ) + 10395*r^2*η^6*sin(θ)^2*cos(θ) - 4620*r^2*η^6*cos(θ) - 18711*r^2*η^4*sin(θ)^2*cos(θ) + 4158*r^2*η^4*cos(θ) + 3465*r^2*η^2*sin(θ)^2*cos(θ) + 1530*r^2*sin(θ)^2*cos(θ) - 452*r^2*cos(θ) + 97020*r*η^4*cos(2*θ) - 41580*r*cos(2*θ) - 1940400*cos(θ))/646800)
end 

```




    Ψ_ηηθ_smallr




```julia
const N = 128;  # Total number of intervals (steps-1) in the η direction
const M = 360;  # Total number of intervals (steps-1) in the θ direction
ψ_j = zeros(N+1,M+1); # Initialise the matrix for the streamfunction at the current time step
const Δη = 1.0/N;
const Δθ = pi/M;

const μ = 1.0
const crank_nic = 1.0 #crank_nic = 1 means implicit and crank_nic =0 means explicit θ-scheme
matrix_condition_numb = zeros(M+1);

meshsize = SA{Int64}[N,M] 
delta_sizes= SA{Float64}[Δθ, Δη] 
paramtrs =  (μ, "3rdOrder")

rey_number = 10.0
```




    10.0




```julia
d3  =  1.0*crank_nic/rey_number        #Diagonal entries for the matrix problem
d2  = -4.0*crank_nic/rey_number 
d1  =  6.0*crank_nic/rey_number

θ_order = "1stOrder"
c    =  Δη/Δθ
h    = 0
matrix_save = 0
for m = 2:M+1 
    rhs = @view ψ_j[2:N,m]
    if m > 3 && θ_order == "2ndOrder"
        c    =  Δη/(2*Δθ)
        h    =  1
    end
    
    A =spdiagm(-2=>zeros(N-3), -1 => zeros(N-2), 0 => zeros(N-1), 1 => zeros(N-2), 2 => zeros(N-3)) 
    row = 1
    n = row + 1
    #A[row,row-2]  = nuemann boundary condition applied and add 0 to rhs
    #A[row,row-1] = 0 dirichlet boundary condition applied, this term is zero
    A[row,row]    = -d3 + d1 + ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*c +2*ψ_η(n, m-1, ψ_j, meshsize)*c
    A[row,row+1]  = d2 -ψ_η(n, m-1, ψ_j, meshsize)*c
    A[row,row+2]  = d3
    
    rhs[row]      = ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*ψ_j[n, m-1-h]*c +
                        -ψ_η(n, m-1, ψ_j, meshsize)*ψ_ηη(n, m-1-h, ψ_j, meshsize)*c +
                    -ψ_ηηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*(1.0-crank_nic)/rey_number
    
    
    row = 2
    n = row + 1
    #A[row,row-2]  = d3 dirichlet boundary condition applied, this term is zero
    A[row,row-1]  = d2 - ψ_η(n, m-1, ψ_j, meshsize)*c 
    A[row,row]    = d1 + ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*c +2*ψ_η(n, m-1, ψ_j, meshsize)*c
    A[row,row+1]  = d2 -ψ_η(n, m-1, ψ_j, meshsize)*c
    A[row,row+2]  = d3
    
    rhs[row]      = ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*ψ_j[n, m-1-h]*c +
                    -ψ_η(n, m-1, ψ_j, meshsize)*ψ_ηη(n, m-1-h, ψ_j, meshsize)*c +
                    -ψ_ηηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*(1.0-crank_nic)/rey_number
                    

    for row = 3:N-3
        n = row + 1
        A[row,row-2]  = d3 + ψ_η(n, m-1, ψ_j, meshsize)*c/12.0
        A[row,row-1]  = d2 - 4*ψ_η(n, m-1, ψ_j, meshsize)*c/3.0 
        A[row,row]    = d1 + ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*c +2.5*ψ_η(n, m-1, ψ_j, meshsize)*c
        A[row,row+1]  = d2 - 4*ψ_η(n, m-1, ψ_j, meshsize)*c/3.0
        A[row,row+2]  = d3 + ψ_η(n, m-1, ψ_j, meshsize)*c/12.0
    
        rhs[row]      = ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*ψ_j[n, m-1-h]*c +
                        -ψ_η(n, m-1, ψ_j, meshsize)*ψ_ηη(n, m-1-h, ψ_j, meshsize)*c+
                        -ψ_ηηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*(1.0-crank_nic)/rey_number
    end

    row = N-2
    n = row + 1
    A[row,row-2]  = d3
    A[row,row-1]  = d2 - ψ_η(n, m-1, ψ_j, meshsize)*c 
    A[row,row]    = d1 + ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*c +2*ψ_η(n, m-1, ψ_j, meshsize)*c
    A[row,row+1]  = d2 -ψ_η(n, m-1, ψ_j, meshsize)*c
    #A[row,row+2]  = d3 dirichlet boundary condition applied, this term is zero
    
    rhs[row]      = ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*ψ_j[n, m-1-h]*c +
                        -ψ_η(n, m-1, ψ_j, meshsize)*ψ_ηη(n, m-1-h, ψ_j, meshsize)*c +
                    -ψ_ηηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*(1.0-crank_nic)/rey_number
    
    row = N-1
    n = row + 1
    A[row,row-2]  = d3
    A[row,row-1]  = d2 - ψ_η(n, m-1, ψ_j, meshsize)*c 
    A[row,row]    = d3 + d1 + ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*c +2*ψ_η(n, m-1, ψ_j, meshsize)*c
    #A[row,row+1]  = d2 -ψ_η(n, m-1, ψ_j, meshsize)*c  dirichlet boundary condition applied, this term is zero
    #A[row,row+2]  = d3 nuemann boundary condition applied and add +2Δημsin(θ)*d3 to rhs
    
    rhs[row]      = ψ_ηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*ψ_j[n, m-1-h]*c +
                        -ψ_η(n, m-1, ψ_j, meshsize)*ψ_ηη(n, m-1-h, ψ_j, meshsize)*c + 
                    -ψ_ηηηη(n, m-1, ψ_j, meshsize, delta_sizes, paramtrs)*(1.0-crank_nic)/rey_number +
                    + 2*Δη*μ*sin(Δθ*(m-1))*d3
    
    

    #LUA = lu(A, check=true)
    #LUA2= factorize(A)
    rhs .= A\rhs
    #ldiv!(LUA2,rhs)
    matrix_condition_numb[m] = cond(Matrix(A),2)
           #This frees up the memory on the next automated garbage collection GC.gc() 
    if m ==37
        matrix_save = A
    end
    A = nothing
end

```


```julia
pyplot()
psi_sol = reshape(ψ_j,length(ψ_j));
η_mesh = [(-1 +(n-1)*Δη)  for m=1:M+1 for n=1:N+1];
θ_mesh = [((m-1)*Δθ) for m=1:M+1 for n=1:N+1 ];
η_vector = range(-1,stop = 1,length = N+1);
θ_vector = range(0,stop = pi,length = M+1);
m_at_fail = M+1
for m = 1:M+1
    if abs(findmax(abs.(ψ_j[:,m]))[1]) > 0.4
        print("Convergence fails at m = $m \ni.e near θ = $(m-1)π/$M\u1D9C radians $(round((m-1)*Δθ,digits=1))\u1D9C \nθ = $(round((m-1)*180/M,digits=1))\U000B0 degrees R_s = $rey_number")
        m_at_fail = m
        break
    elseif m == M+1
        println("Solution seems to be smooth")
    end
    
end
```

    Convergence fails at m = 99 
    i.e near θ = 98π/360ᶜ radians 0.9ᶜ 
    θ = 49.0° degrees R_s = 10.0


```julia
compare_with_timemarch = true
ψ_tm = 0;
if compare_with_timemarch# Load file with matrix N by M 
    try 
    ψ_tm =readdlm("psi_solution_R_45.0_1stOrder_sym.txt",Float64,); #Check that N and M dimensions match
    if size(ψ_tm) == (N+1,M+1) 
        print("Meshsize dimensions of the timemarched ψ match with the θ-marched solution \nthis is good!")
    else
        print("Dimensions of timemarched ψ DO NOT match with the θ-marched  which are $(size(ψ_tm))\n ERROR..")
        compare_with_timemarch = false;
    end
    catch e
        println("Could not open file don't worry code will continue regardless")
        compare_with_timemarch = false;
    end
end
```

    Meshsize dimensions of the timemarched ψ match with the θ-marched solution 
    this is good!


```julia

wireframe( range(-1.0,stop = 1.0,length = round(Int,N/4)+1),range(0.0,stop = pi,length = round(Int,M/4)+1), transpose(ψ_j[1:4:end,1:4:end]), colour=:greys, fmt = :svg)
yticks!(([0:π/6:π;], ["0", "\\pi/6", "\\pi/3", "\\pi/2", "2\\pi/3", "5\\pi/6","\\pi"]))
plot!(xlabel=L"\eta", ylabel=L"\theta", zlabel=L"\Psi", title=L"R_s = %$rey_number"*"\n"*L"N_η=%$N, M_θ=%$M", size=(1000,800),camera=(-60,30) )
plot!(xtickfont = font(12), ytickfont = font(12), ztickfont = font(12), titlefont = font(11))
#yticks!(([0:π/6:π;], ["0", "\\pi/6", "\\pi/3", "\\pi/2", "2\\pi/3", "5\\pi/6","\\pi"]))
#plot!(xlabel=L"η", ylabel=L"θ", zlabel=L"Ψ", title=L"R_s = %$rey_number"*"\n"*L"N_η=%$N, M_θ=%$M")
#plot!(xtickfont = font(12), ytickfont = font(12), ztickfont = font(12), titlefont=font(10))
```




    

    




```julia

m=m_at_fail - 7
m = 90

eta=range(0, stop=1,length= N+1)
theta_march_plt = plot(eta, ψ_j[:,m], label="θ-marched numerical \n steady solution", linecolor=:red,    linestyle=:dot,   linewidth =2, legend=:outertopright, size=(1100,500))
xlabel!("η \n"*L"N_η=%$N, M_θ=%$M")
ylabel!(L"Ψ", yguidefontrotation=-90)
title!("η-profile  "*L"θ=\frac{%$(m-1)}{%$M}π \quad R_s = %$rey_number")
if compare_with_timemarch
    plot!(eta, ψ_tm[:,m], label="timemarched numerical \nsolution", linecolor=:cyan,    linestyle=:dot,   linewidth =2, legend=:outertopright, size=(1100,500))
end
if rey_number < 60
 plot!(x->Ψ_smallr(Δθ*(m-1),x,rey_number,3),0,1, label="Small "*L"R_s"*" asymptotic",linestyle=:dot,   linewidth =3, linecolor=:blue)
end

display(theta_march_plt)

ψ_η_plt = plot(eta[2:N], [ψ_η(n, m, ψ_j, meshsize)/Δη for n=2:N ], title=L"ψ_η",label="θ-marched" , linestyle=:dot, linewidth=2, linecolor=:red, size=(1100,500))
ψ_ηη_plt = plot(eta[2:N], [ψ_ηη(n, m, ψ_j, meshsize)/Δη^2 for n=2:N ], title=L"ψ_{ηη}",label="θ-marched" , linestyle=:dot, linewidth=2,linecolor=:red,  size=(1100,500))
ψ_ηηη_plt = plot(eta[2:N], [ψ_ηηη(n, m, ψ_j, meshsize,  delta_sizes, paramtrs)/Δη^3 for n=2:N ], title=L"ψ_{ηηη}",linewidth=2,linecolor=:red, label="θ-marched" , linestyle=:dot, size=(1100,500))
ψ_ηηηη_plt = plot(eta[2:N], [ψ_ηηηη(n, m, ψ_j, meshsize,  delta_sizes, paramtrs)/Δη^4 for n=2:N ], title=L"ψ_{ηηηη}",linewidth=2,linecolor=:red, label="θ-marched" , linestyle=:dot, size=(1100,500))
ψ_θ_plt = plot(eta[2:N], [(ψ_j[n,m+1]-ψ_j[n,m])/Δθ for n=2:N ], title=L"ψ_{θ}",label="θ-marched" , linestyle=:dot, linewidth=2, linecolor=:red,size=(1100,500))
ψ_ηηθ_plt = plot(eta[2:N], [(ψ_ηη(n, m+1, ψ_j, meshsize)-ψ_ηη(n, m, ψ_j, meshsize))/(Δθ*Δη^2) for n=2:N ], title=L"ψ_{ηηθ}",label="θ-marched" , linecolor=:red,linewidth=2, linestyle=:dot, size=(1100,500))
if rey_number < 50.0
    ψ_η_plt =  plot(ψ_η_plt,    x-> Ψ_η_smallr(Δθ*(m-1), x ,rey_number), 0.0,1.0, linecolor=:blue, label="Asymptotic", legend=:outertopright)
   ψ_ηη_plt =  plot(ψ_ηη_plt,   x-> Ψ_ηη_smallr(Δθ*(m-1), x ,rey_number), 0.0,1.0, linecolor=:blue,label="Asymptotic",legend=false)
  ψ_ηηη_plt =  plot(ψ_ηηη_plt,  x-> Ψ_ηηη_smallr(Δθ*(m-1), x ,rey_number), 0.0,1.0, linecolor=:blue,label="Asymptotic", legend=:outertopright)
 ψ_ηηηη_plt =  plot(ψ_ηηηη_plt, x-> Ψ_ηηηη_smallr(Δθ*(m-1), x ,rey_number), 0.0,1.0, linecolor=:blue,label="Asymptotic", legend=false)
  ψ_ηηθ_plt =  plot(ψ_ηηθ_plt,   x-> Ψ_ηηθ_smallr(Δθ*(m-1), x ,rey_number), 0.0,1.0, linecolor=:blue, label="Asymptotic", legend=false)
  ψ_θ_plt =  plot(ψ_θ_plt,   x-> Ψ_θ_smallr(Δθ*(m-1), x ,rey_number), 0.0,1.0, linecolor=:blue, label="Asymptotic",legend=:outertopright)  
end

plot(ψ_η_plt, ψ_ηη_plt, ψ_ηηη_plt, ψ_ηηηη_plt, ψ_θ_plt, ψ_ηηθ_plt, layout = (3, 2), size=(1100,1000))
xlabel!("η at θ=$(round(Δθ*(m-1),digits=1))\u1D9C")
ylabel!(L"Ψ", yguidefontrotation=-90)







```


```julia
eta_profile =0.6
theta=range(0, stop=pi,length= M+1)

plot((0:m_at_fail-1)*pi/(M+1),ψ_j[round(Int,(eta_profile)/Δη)+1,1:m_at_fail], linecolor=:red, label="θ-marched numerical \n steady solution", linestyle=:dot, linewidth=1,legend=:outertopright, size=(1100,500))
if rey_number < 60
plot!(x->Ψ_smallr(x,eta_profile,rey_number, 3),0,pi, label="Small "*L"R_s"*" asymptotic", linestyle=:dot, linecolor=:blue, linewidth=3)
end
xticks!(([0:π/6:π;], ["0", "\\pi/6", "\\pi/3", "\\pi/2", "2\\pi/3", "5\\pi/6","\\pi"]))
title!(L"θ"*"-profile at "*L"η= %$eta_profile, \ R_s = %$rey_number")
xlabel!("θ \n"*L"N_η=%$N, M_θ=%$M")
ylabel!(L"Ψ", yguidefontrotation=-90)
```


```julia
plot(theta,log.(matrix_condition_numb), label="Log of Matrix \ncondition number",  size=(1100,500), legend=:outertopright)
ylabel!("log(cond(A)")
xlabel!("θ")
xticks!(([0:π/6:π;], ["0", "\\pi/6", "\\pi/3", "\\pi/2", "2\\pi/3", "5\\pi/6","\\pi"]))
```


```julia
t = (1.0,"Order")
```




    (1.0, "Order")




```julia
p=1
println(matrix_save[p,p:p+3]')
p=2
println(matrix_save[p,p-1:p+3]')
for p =3:N-3
println(matrix_save[p,p-2:p+2]')
end
p=N-2
println(matrix_save[p,p-2:p+1]')
p=N-1
println(matrix_save[p,p-2:p]')
```

    [0.5023953221103209 -0.4011978953921838 0.1 0.0]
    [-0.40119727049317006 0.6023940723131727 -0.40119727049317006 0.1 0.0]
    [0.10009967490070991 -0.4015947984113585 0.6029897783495145 -0.4015947984113585 0.10009967490070991]
    [0.10009953820434786 -0.4015926112695657 0.6029856774612735 -0.4015926112695657 0.10009953820434786]
    [0.10009936245223258 -0.4015897992357213 0.6029804049018875 -0.4015897992357213 0.10009936245223258]
    [0.10009914764470709 -0.40158636231531336 0.6029739606821726 -0.40158636231531336 0.10009914764470709]
    [0.10009889378227999 -0.40158230051647975 0.6029663448180285 -0.40158230051647975 0.10009889378227999]
    [0.10009860086567894 -0.4015776138508629 0.6029575573320405 -0.4015776138508629 0.10009860086567894]
    [0.10009826889591361 -0.40157230233461766 0.6029475982553687 -0.40157230233461766 0.10009826889591361]
    [0.10009789787434815 -0.40156636598957046 0.6029364676299195 -0.40156636598957046 0.10009789787434815]
    [0.10009748780278303 -0.4015598048445284 0.6029241655107992 -0.4015598048445284 0.10009748780278303]
    [0.10009703868354602 -0.4015526189367364 0.6029106919690443 -0.4015526189367364 0.10009703868354602]
    [0.10009655051959253 -0.4015448083134803 0.602896047094627 -0.4015448083134803 0.10009655051959253]
    [0.10009602331461467 -0.40153637303383466 0.602880230999729 -0.40153637303383466 0.10009602331461467]
    [0.10009545707315949 -0.4015273131705519 0.6028632438222827 -0.4015273131705519 0.10009545707315949]
    [0.10009485180075568 -0.40151762881209097 0.6028450857297729 -0.40151762881209097 0.10009485180075568]
    [0.10009420750404893 -0.40150732006478274 0.602825756923295 -0.40150732006478274 0.10009420750404893]
    [0.10009352419094561 -0.40149638705512974 0.6028052576418641 -0.40149638705512974 0.10009352419094561]
    [0.10009280187076477 -0.40148482993223633 0.6027835881669702 -0.40148482993223633 0.10009280187076477]
    [0.10009204055439799 -0.40147264887036777 0.6027607488273734 -0.40147264887036777 0.10009204055439799]
    [0.1000912402544771 -0.40145984407163354 0.6027367400041332 -0.40145984407163354 0.1000912402544771]
    [0.10009040098554954 -0.40144641576879253 0.6027115621358652 -0.40144641576879253 0.10009040098554954]
    [0.10008952276426103 -0.4014323642281763 0.6026852157242197 -0.4014323642281763 0.10008952276426103]
    [0.10008860560954545 -0.4014176897527271 0.6026577013395755 -0.4014176897527271 0.10008860560954545]
    [0.10008764954282169 -0.40140239268514694 0.602629019626941 -0.40140239268514694 0.10008764954282169]
    [0.10008665458819711 -0.40138647341115363 0.6025991713120578 -0.40138647341115363 0.10008665458819711]
    [0.10008562077267757 -0.40136993236284113 0.6025681572076969 -0.40136993236284113 0.10008562077267757]
    [0.10008454812638369 -0.4013527700221389 0.6025359782201429 -0.4013527700221389 0.10008454812638369]
    [0.100083436682773 -0.40133498692436786 0.6025026353558584 -0.40133498692436786 0.100083436682773]
    [0.10008228647886797 -0.4013165836618874 0.6024681297283188 -0.4013165836618874 0.10008228647886797]
    [0.10008109755548943 -0.4012975608878308 0.6024324625650143 -0.4012975608878308 0.10008109755548943]
    [0.10007986995749532 -0.40127791931992507 0.602395635214607 -0.40127791931992507 0.10007986995749532]
    [0.10007860373402434 -0.40125765974438937 0.60235764915424 -0.40125765974438937 0.10007860373402434]
    [0.10007729893874442 -0.40123678301991067 0.602318505996986 -0.40123678301991067 0.10007729893874442]
    [0.10007595563010563 -0.40121529008169005 0.6022782074994341 -0.40121529008169005 0.10007595563010563]
    [0.10007457387159738 -0.40119318194555786 0.6022367555693996 -0.40119318194555786 0.10007457387159738]
    [0.10007315373200949 -0.4011704597121519 0.6021941522737578 -0.4011704597121519 0.10007315373200949]
    [0.10007169528569723 -0.40114712457115576 0.602150399846386 -0.40114712457115576 0.10007169528569723]
    [0.10007019861284956 -0.40112317780559287 0.6021055006962132 -0.40112317780559287 0.10007019861284956]
    [0.1000686637997608 -0.4010986207961728 0.6020594574153655 -0.4010986207961728 0.1000686637997608]
    [0.10006709093910532 -0.4010734550256851 0.6020122727874022 -0.4010734550256851 0.10006709093910532]
    [0.10006548013021488 -0.40104768208343794 0.6019639497956358 -0.40104768208343794 0.10006548013021488]
    [0.10006383147935857 -0.4010213036697371 0.6019144916315271 -0.4010213036697371 0.10006383147935857]
    [0.10006214510002513 -0.4009943216004021 0.6018639017031495 -0.4009943216004021 0.10006214510002513]
    [0.10006042111320719 -0.400966737811315 0.6018121836437161 -0.400966737811315 0.10006042111320719]
    [0.10005865964768748 -0.4009385543629996 0.6017593413201637 -0.4009385543629996 0.10005865964768748]
    [0.10005686084032668 -0.4009097734452269 0.6017053788417844 -0.4009097734452269 0.10005686084032668]
    [0.10005502483635273 -0.4008803973816436 0.6016503005689025 -0.4008803973816436 0.10005502483635273]
    [0.1000531517896513 -0.40085042863442066 0.6015941111215869 -0.40085042863442066 0.1000531517896513]
    [0.10005124186305737 -0.4008198698089179 0.6015368153883971 -0.4008198698089179 0.10005124186305737]
    [0.1000492952286477 -0.4007887236583631 0.601478418535152 -0.4007887236583631 0.1000492952286477]
    [0.10004731206803384 -0.40075699308854135 0.6014189260137225 -0.40075699308854135 0.10004731206803384]
    [0.1000452925726558 -0.4007246811624929 0.6013583435708363 -0.4007246811624929 0.1000452925726558]
    [0.10004323694407599 -0.4006917911052159 0.6012966772568943 -0.4006917911052159 0.10004323694407599]
    [0.10004114539427328 -0.4006583263083724 0.6012339334347926 -0.4006583263083724 0.10004114539427328]
    [0.10003901814593714 -0.4006242903349942 0.6011701187887443 -0.4006242903349942 0.10003901814593714]
    [0.10003685543276171 -0.4005896869241874 0.6011052403330982 -0.4005896869241874 0.10003685543276171]
    [0.1000346574997395 -0.40055451999583197 0.60103930542115 -0.40055451999583197 0.1000346574997395]
    [0.10003242460345477 -0.40051879365527615 0.600972321753942 -0.40051879365527615 0.10003242460345477]
    [0.10003015701237636 -0.4004825121980216 0.6009042973890467 -0.4004825121980216 0.10003015701237636]
    [0.10002785500714999 -0.4004456801143997 0.6008352407493337 -0.4004456801143997 0.10002785500714999]
    [0.10002551888088972 -0.40040830209423534 0.600765160631713 -0.40040830209423534 0.10002551888088972]
    [0.10002314893946865 -0.40037038303149824 0.6006940662158553 -0.40037038303149824 0.10002314893946865]
    [0.10002074550180871 -0.4003319280289394 0.6006219670728843 -0.4003319280289394 0.10002074550180871]
    [0.10001830890016947 -0.40029294240271146 0.6005488731740389 -0.40029294240271146 0.10001830890016947]
    [0.10001583948043578 -0.40025343168697247 0.6004747948993044 -0.40025343168697247 0.10001583948043578]
    [0.10001333760240445 -0.4002134016384711 0.6003997430460094 -0.4002134016384711 0.10001333760240445]
    [0.10001080364006956 -0.40017285824111287 0.6003237288373867 -0.40017285824111287 0.10001080364006956]
    [0.1000082379819067 -0.40013180771050705 0.6002467639310983 -0.40013180771050705 0.1000082379819067]
    [0.10000564103115578 -0.40009025649849234 0.6001688604277202 -0.40009025649849234 0.10000564103115578]
    [0.10000301320610261 -0.40004821129764156 0.6000900308791891 -0.40004821129764156 0.10000301320610261]
    [0.10000035494035908 -0.40000567904574535 0.6000102882972083 -0.40000567904574535 0.10000035494035908]
    [0.09999766668314207 -0.399962666930273 0.5999296461616122 -0.399962666930273 0.09999766668314207]
    [0.09999494889955075 -0.39991918239281204 0.5998481184286902 -0.39991918239281204 0.09999494889955075]
    [0.09999220207084282 -0.39987523313348505 0.5997657195394679 -0.39987523313348505 0.09999220207084282]
    [0.09998942669470902 -0.39983082711534423 0.5996824644279478 -0.39983082711534423 0.09998942669470902]
    [0.09998662328554651 -0.39978597256874415 0.5995983685293078 -0.39978597256874415 0.09998662328554651]
    [0.09998379237473078 -0.3997406779956924 0.5995134477880585 -0.3997406779956924 0.09998379237473078]
    [0.09998093451088615 -0.3996949521741783 0.5994277186661607 -0.3996949521741783 0.09998093451088615]
    [0.09997805026015504 -0.39964880416248055 0.5993411981511023 -0.39964880416248055 0.09997805026015504]
    [0.09997514020646592 -0.3996022433034545 0.599253903763937 -0.3996022433034545 0.09997514020646592]
    [0.09997220495179993 -0.3995552792287988 0.5991658535672854 -0.3995552792287988 0.09997220495179993]
    [0.09996924511645644 -0.3995079218633029 0.5990770661732991 -0.3995079218633029 0.09996924511645644]
    [0.09996626133931726 -0.399460181429076 0.5989875607515917 -0.399460181429076 0.09996626133931726]
    [0.0999632542781099 -0.39941206844975824 0.5988973570371335 -0.39941206844975824 0.0999632542781099]
    [0.09996022460966968 -0.39936359375471475 0.5988064753381189 -0.39936359375471475 0.09996022460966968]
    [0.09995717303020096 -0.39931476848321534 0.598714936543802 -0.39931476848321534 0.09995717303020096]
    [0.09995410025553741 -0.39926560408859857 0.598622762132306 -0.39926560408859857 0.09995410025553741]
    [0.09995100702140146 -0.3992161123424233 0.5985299741784086 -0.3992161123424233 0.09995100702140146]
    [0.09994789408366304 -0.3991663053386086 0.5984365953613048 -0.3991663053386086 0.09994789408366304]
    [0.09994476221859769 -0.399116195497563 0.5983426489723517 -0.399116195497563 0.09994476221859769]
    [0.09994161222314411 -0.39906579557030586 0.5982481589227957 -0.39906579557030586 0.09994161222314411]
    [0.09993844491516132 -0.3990151186425812 0.598153149751489 -0.3990151186425812 0.09993844491516132]
    [0.09993526113368546 -0.3989641781389673 0.5980576466325948 -0.3989641781389673 0.09993526113368546]
    [0.09993206173918638 -0.3989129878269819 0.5979616753832877 -0.3989129878269819 0.09993206173918638]
    [0.09992884761382426 -0.3988615618211881 0.5978652624714527 -0.3988615618211881 0.09992884761382426]
    [0.09992561966170622 -0.3988099145872996 0.5977684350233854 -0.3988099145872996 0.09992561966170622]
    [0.09992237880914318 -0.3987580609462908 0.5976712208314989 -0.3987580609462908 0.09992237880914318]
    [0.09991912600490702 -0.39870601607851214 0.5975736483620431 -0.39870601607851214 0.09991912600490702]
    [0.09991586222048837 -0.39865379552781377 0.5974757467628384 -0.39865379552781377 0.09991586222048837]
    [0.09991258845035508 -0.39860141520568126 0.5973775458710323 -0.39860141520568126 0.09991258845035508]
    [0.09990930571221156 -0.39854889139538496 0.5972790762208814 -0.39854889139538496 0.09990930571221156]
    [0.09990601504725918 -0.3984962407561469 0.597180369051568 -0.3984962407561469 0.09990601504725918]
    [0.09990271752045805 -0.39844348032732874 0.5970814563150534 -0.39844348032732874 0.09990271752045805]
    [0.09989941422079018 -0.398390627532643 0.5969823706839782 -0.398390627532643 0.09989941422079018]
    [0.09989610626152458 -0.39833770018439324 0.5968831455596154 -0.39833770018439324 0.09989610626152458]
    [0.09989279478048413 -0.398284716487746 0.596783815079883 -0.398284716487746 0.09989279478048413]
    [0.0998894809403149 -0.39823169504503836 0.596684414127424 -0.39823169504503836 0.0998894809403149]
    [0.09988616592875794 -0.3981786548601271 0.5965849783377657 -0.3981786548601271 0.09988616592875794]
    [0.09988285095892395 -0.39812561534278323 0.5964855441075622 -0.39812561534278323 0.09988285095892395]
    [0.09987953726957106 -0.39807259631313696 0.5963861486029339 -0.39807259631313696 0.09987953726957106]
    [0.09987622612538626 -0.3980196180061801 0.5962868297679131 -0.3980196180061801 0.09987622612538626]
    [0.09987291881727062 -0.39796670107632987 0.5961876263330066 -0.39796670107632987 0.09987291881727062]
    [0.0998696166626289 -0.3979138666020623 0.5960885778238881 -0.3979138666020623 0.0998696166626289]
    [0.09986632100566382 -0.397861136090621 0.5959897245702326 -0.397861136090621 0.09986632100566382]
    [0.09986303321767556 -0.397808531482809 0.5958911077147075 -0.397808531482809 0.09986303321767556]
    [0.09985975469736698 -0.3977560751578715 0.5957927692221341 -0.3977560751578715 0.09985975469736698]
    [0.09985648687115492 -0.39770378993847866 0.5956947518888364 -0.39770378993847866 0.09985648687115492]
    [0.09985323119348855 -0.39765169909581666 0.5955970993521952 -0.39765169909581666 0.09985323119348855]
    [0.09984998914717484 -0.3975998263547975 0.5954998561004244 -0.3975998263547975 0.09984998914717484]
    [0.09984676224371239 -0.3975481958993982 0.5954030674825901 -0.3975481958993982 0.09984676224371239]
    [0.0998435520236338 -0.3974968323781407 0.5953067797188946 -0.3974968323781407 0.0998435520236338]
    [0.09984036005685781 -0.39744576090972494 0.5952110399112491 -0.39744576090972494 0.09984036005685781]
    [0.09983718794305178 -0.3973950070888284 0.5951158960541596 -0.3973950070888284 0.09983718794305178]
    [0.09983403731200548 -0.3973445969920877 0.5950213970501588 -0.3973445969920877 0.09983403731200548]
    [0.1 -0.39797091789101047 0.5959421338389768 -0.39797091789101047]
    [0.1 -0.39793373915081115 0.6958677970121524]



```julia
det(matrix_save)
```




    9.445352039164064e-121



### Finte difference stencil generator
See https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_finite_difference 


```julia
s = [-3, -2, -1, 0, 1] #Stencil points at which finite differences will be calculated
derivs_required = 1
s_N = length(s)
a_vector = zeros(s_N)

stencil_matrix = zeros(s_N,s_N)
for i=1:s_N
    for j=1:s_N
        stencil_matrix[i,j] = s[j]^(i-1)
    end
end

a_vector[derivs_required+1] = 1*factorial(derivs_required)

(stencil_matrix\a_vector)'
```




    1×5 adjoint(::Vector{Float64}) with eltype Float64:
     -0.0833333  0.5  -1.5  0.833333  0.25




```julia


```


```julia
factorial(3)
```




    6




```julia

```

b8f19c78b77680df49ba5c0c263333a1e144807e
ipynb
Jupyter Notebook
Psi_SteadyStates_ThetaMarching_SYM.ipynb
gsagoo/SolvingDifferentialEquations
63c2b231b45ad64206a11824198b3d5fe20c062b
Unlicense
Psi_SteadyStates_ThetaMarching_SYM.ipynb
gsagoo/SolvingDifferentialEquations
63c2b231b45ad64206a11824198b3d5fe20c062b
Unlicense
Psi_SteadyStates_ThetaMarching_SYM.ipynb
gsagoo/SolvingDifferentialEquations
63c2b231b45ad64206a11824198b3d5fe20c062b
Unlicense
Qwen/Qwen-72B
1. YES
2. YES