# Did i coded my pde correctly?

**URL:** https://discourse.julialang.org/t/did-i-coded-my-pde-correctly/121661
**Category:** Numerics
**Tags:** pde
**Created:** [October 23, 2024, 5:08pm UTC](https://discourse.julialang.org/t/did-i-coded-my-pde-correctly/121661 "2024-10-23T17:08:55Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![Shashi2\_Kumar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shashi2_kumar/32/212826_2.png) [@Shashi2\_Kumar](https://discourse.julialang.org/u/Shashi2_Kumar)
#### Post date: [October 23, 2024, 5:08pm UTC](https://discourse.julialang.org/t/did-i-coded-my-pde-correctly/121661/1 "2024-10-23T17:08:55Z")

</div>

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?

 ![EC](https://global.discourse-cdn.com/julialang/original/3X/4/9/49508216de77b1974a13d272d6e838dab9df0480.png)

# 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)