After many weeks of search for a pde solver for my elliptic pde i came across Julia. Found that it has a FNC package (FundamentalsNumericalComputation) that solves pde. I tried to encode my pde with boundary conditions as in image, but not sure if i did it correctly. Especially have doubt about boundary conditions. I am learning and solving at the same time. Can you guys tell mistake?

# Parameters

n = 3 # polytropic index

m = 40 # points in ξ direction (radial)

k = 40 # points in μ direction (angular)

xiE = 6.897 # for n=3 polytrope

v = 0 # Dimensionless rotation parameter (as given)

# Define the PDE operator

function ϕ(X, Y, U, Ux, Uxx, Uy, Uyy)

ξ = X # X represents ξ grid

μ = Y # Y represents μ grid

```
# Handle the boundary conditions
ξ_boundary = @. (X ≈ 0)
μ_boundary_lower = @. (Y ≈ -1)
μ_boundary_upper = @. (Y ≈ 1)
# First term: radial part
term1 = @. (1 / ξ^2) * (2 * ξ * Ux + ξ^2 * Uxx) # Expanded form of ∂ξ(ξ² ∂ξ(θ))
# Second term: angular part
term2 = @. (1 / ξ^2) * ((1 - μ^2) * Uyy - 2 * μ * Uy) # Expanded form of ∂μ((1-μ²) ∂μ(θ))
# The right-hand side with the polytropic index
R = @. term1 + term2 - v + U^n
# Apply boundary conditions at ξ = 0
# Dirichlet condition: θ(0, μ) = 1
R[ξ_boundary] = @. U[ξ_boundary] - 1
# Neumann condition: ∂θ/∂ξ (0, μ) = 0 at ξ = 0
# Already enforced at ξ = 0 due to the equation structure
R[ξ_boundary] = @. Ux[ξ_boundary]
# Neumann conditions at μ = ±1 (poles)
R[μ_boundary_lower] = @. Uy[μ_boundary_lower]
R[μ_boundary_upper] = @. Uy[μ_boundary_upper]
return R
```

end

# Boundary conditions - Dirichlet at ξ = 0

g(x,y) = (x ≈ 0) ? 1.0 : 0.0 # θ(0,μ) = 1

# Solve using the provided elliptic solver

u = FNC.elliptic(ϕ, g, m, [0,xiE], k, [-1,1])

# Plot the solution

using Plots

# Create a grid for plotting

ξ = range(0, xiE, length=100)

μ = range(-1, 1, length=100)

U = [u(x,y) for x in ξ, y in μ]

# Plot the solution

contourf(ξ, μ, U’,

xlabel=“ξ”,

ylabel=“μ”,

title=“Emden-Chandrasekhar Solution (n=$n)”,

color=:viridis,

aspect_ratio=:equal)