# Problem after discretizing a PDE to ODE and solve with DifferentialEquations.jl

**URL:** <https://discourse.julialang.org/t/problem-after-discretizing-a-pde-to-ode-and-solve-with-differentialequations-jl/81451>\
**Category:** General Usage\
**Tags:** differentialequation\
**Created:** [May 22, 2022, 6:53am UTC](https://discourse.julialang.org/t/problem-after-discretizing-a-pde-to-ode-and-solve-with-differentialequations-jl/81451 "2022-05-22T06:53:24Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Molox](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/molox/32/35002_2.png) [@Molox](https://discourse.julialang.org/u/Molox)\
**Post date:** [May 22, 2022, 6:53am UTC](https://discourse.julialang.org/t/problem-after-discretizing-a-pde-to-ode-and-solve-with-differentialequations-jl/81451/1 "2022-05-22T06:53:24Z")

</div>

Hi, I was reading [https://www.stochasticlifestyle.com/solving-systems-stochastic-pdes-using-gpus-julia/](https://www.stochasticlifestyle.com/solving-systems-stochastic-pdes-using-gpus-julia/)

And I would like to know if it is possible to discretize a 4th derivative PDE to ODE: ddt(u) = ddddx(u) and solve with DifferentialEquations.jl.  
I am trying to discretize the RHS of the equation and solve it. But I got an error with **Instability detected. Aborting.** I increased the mesh size and changed to different solvers, but it is still not working.

Any help would be appreciated!

using DiffEqOperators  
using BandedMatrices  
using DifferentialEquations

m = 0.5  
L = pi/m  
n = 200 # 0,1,…,n-1,n  
N = n-2 #2,3,…,n-3,n-2  
# mesh from 0 to n  
xl = LinRange(0.0,L,n)  
dx = xl[2] - xl[1]

A = BandedMatrix(  
-2 =\> ones(N-2)_(1),  
-1 =\> ones(N-1)_(-4),  
0 =\> ones(N)_(6),  
1 =\> ones(N-1)_(-4),  
2 =\> ones(N-2)\*(1)  
)

A[1,1] = 1  
A[1,2] = -4/5  
A[1,3] = 1/5  
A[end,end] = 1  
A[end,end-1] = -4/5  
A[end,end-2] = 1/5

function myode(ddu, du, u, p, t)  
ddu[:] = A\*u./(dx^4)  
end

x0 = LinRange(0, L, n)  
u0 = sin.(m .\*x0)[2:end-1]  
du0 = cos.(m .\*x0)[2:end-1]  
tmax = 1.  
p = 0  
tspan = (0,tmax)  
prob = SecondOrderODEProblem(myode,du0,u0,tspan,p)  
sol = solve(prob, RadauIIA5(), saveat=0:0.01:tmax)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 22, 2022, 9:42am UTC](https://discourse.julialang.org/t/problem-after-discretizing-a-pde-to-ode-and-solve-with-differentialequations-jl/81451/2 "2022-05-22T09:42:01Z")

</div>

Are you sure that’s a stable discretization of the PDE? Do a Von Neumann analysis.

For reference, look through the following:

> [@PSA: How to help yourself debug differential equation solving issues](https://discourse.julialang.org/t/psa-how-to-help-yourself-debug-differential-equation-solving-issues/62489):
>
> Debugging differential equation solver issues almost always boils down to doing the same thing, so this is a summary to help you out. For more information on specific issues, check out the FAQ section of the DifferentialEquations.jl documentation which highlights common issues and questions: [https://diffeq.sciml.ai/dev/basics/faq/](https://diffeq.sciml.ai/dev/basics/faq/)How do I debug why the differential equation solver is diverging? dt \<= dtmin. Aborting. There is either an error in your model specification or the true solution i…