# Implement of Euler Method of Parabolic PDE

I'm experimenting with different schemes for numerically solving PDEs in Mathematica. I use NestList to implement the Euler Method for 1-D PDE, the code is concise:

forwardEuler[t0_, y0_, tend_, h_, func_] :=
Module[{steps = Floor[(tend - t0)/h], ics = {t0, y0}},
NestList[{#[[1]] + h, #[[2]] + h *func /. {t -> #[[1]], y -> #[[2]]}} &,
ics, steps]]


and it works well.

ListLinePlot[{Table[Evaluate[{t,DSolveValue[{y'[t] == (y[t] + t^2 - 2)/(t + 1), y[0] == 2}, y[t],t]}], {t, 0, 6, 0.2}],
forwardEuler[0, 2, 6, 0.2, (y + t^2 - 2)/(t + 1)]},PlotStyle -> {Blue, {Dashed, Orange}}, PlotLegends -> {"Exact Solution", "Forward Euler"}]


However, now I want to use Euler Method and central finite difference to solve PDE like the following.

$$\left\{\begin{array}{l} \frac{\partial u}{\partial t}-\frac{\partial^{2} u}{\partial x^{2}}=x \mathrm{e}^{t}-6 x, \quad 0

Is it still recommended to write an Euler method using NestList (or any function from the Nest family)? It has a simple form in 1-D, but I'd like to see it in 2-D or even 3-D.

Additional: I also write codes of Implicit Euler, Trapezoidal, and Heun for 1D equations:

backwardEuler[t0_, y0_, tend_, h_, func_, tol_] := Module[{func1},
func1[{\[FormalT]_, \[FormalY]_}] := {\[FormalT] + h, NestWhile[N[\[FormalY] + h*func /. {t -> \[FormalT] + h, y -> #}] &, \[FormalY] + h*func /. {t -> \[FormalT], y -> \[FormalY]}, (Abs[#1 - #2] > tol) &, 2, 10]};
NestList[func1[#] &, {t0, y0}, Floor[(tend - t0)/h]]]

trapezoidal[t0_,y0_,tend_,h_,func_,tol_]:=Module[{func1},
func1[{\[FormalT]_,\[FormalY]_}]:={\[FormalT]+h,NestWhile[N[\[FormalY]+h*((func/.{t->\[FormalT],y->\[FormalY]})+(func/.{t->\[FormalT]+h,y->#}))/2]&,\[FormalY]+h*func/.{t->\[FormalT],y->\[FormalY]},(Abs[#1-#2]>tol)&,2,10]};
NestList[func1[#]&,{t0,y0},Floor[(tend-t0)/h]]]

HuenEuler[t0_,y0_,tend_,h_,func_]:=Module[{func1},
func1[{\[FormalT]_,\[FormalY]_}]:={\[FormalT]+h,N[\[FormalY]+h*((func/.{t->\[FormalT],y->\[FormalY]})+(func/.{t->\[FormalT]+h,y->(\[FormalY]+h*func/.{t->\[FormalT],y->\[FormalY]})}))/2]};
NestList[func1[#]&,{t0,y0},Floor[(tend-t0)/h]]]

• Euler method etc. are methods that can be extended to work on ODE system naturally, but your implementation can not. What you really need to think about is: how to code Euler method etc. so they work on ODE system? ( A single ODE is just a special case of ODE system. ) Related: mathematica.stackexchange.com/a/274162/1871 Jan 10 at 10:00