To start with, I need to solve this partial equation numerically, but I do not know how to do that. If I try a finite difference method, I face a problem that $u_{i,j}$ is also located in exponential, consequently, there is no way to find $u_{i,j}$. I tried this fragmentation: $$\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2} + \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{\Delta y^2} = \mathrm{e}^{u_{i,j}}$$ What can help me to solve this equation? What ways do exist to solve this problems?

Boundary condiditions: $$u(0,y)=a, \quad u_x(0,y) = b$$ $$u(x,0)=c, \quad u_y(x,0) = d$$

  • 2
    $\begingroup$ For a beginner-level approach, one could add the time-derivative term there and solve it as a time-evolution problem to the steady state, using the method of lines and a standard ODE integrator. $\endgroup$ Nov 22, 2023 at 23:34
  • 4
    $\begingroup$ You just happen to get a nonlinear system of equations, as shown in your first formula. You solve that like any other nonlinear system, for example with Newton's method. $\endgroup$ Nov 23, 2023 at 0:35
  • 1
    $\begingroup$ If you want more insight into the exact methods, solvers, and preconditioners that are used on this problem, it is called the Bratu problem and is a decent amount of literature on its numerical solution. $\endgroup$
    – whpowell96
    Nov 23, 2023 at 23:06
  • $\begingroup$ You are using two boundary conditions on $x=0$ and two conditions on $y=0$, and nothing on the other parts of the boundary. Is this what you want to do ? $\endgroup$
    – cfdlab
    Nov 28, 2023 at 13:28

1 Answer 1


You can use something like MethodOfLines.jl to solve your problem. You can modify the following example, something like this

using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve

@parameters x y
@variables u(..)

Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ exp(u(x,y))

a = 5.2
b = 4.3
c = 3.5
d = 2.2

bcs = [u(0, y) ~ a,
       u(x, 0) ~ c,
       Dx(u(0, y)) ~ b,
       Dy(u(x, 0)) ~ d]

domains = [x ∈ Interval(0.0, 1.0),
           y ∈ Interval(0.0, 1.0)]

@named pdesys = PDESystem([eq], bcs, domains, [x, y], [u(x, y)])

dx = 0.1
dy = 0.1

discretization = MOLFiniteDifference([x => dx, y => dy], nothing, approx_order=2)

prob = discretize(pdesys, discretization)
sol = NonlinearSolve.solve(prob, NewtonRaphson())

u_sol = sol[u(x, y)]

using Plots

heatmap(sol[x], sol[y], u_sol, xlabel="x", ylabel="y",

Where I've plugged in some random boundary conditions. You might want to tweak some parameters like tolerance, maximum iterations to better suit your problem.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.