I am pleased to announce the official release of SkeelBerzins.jl v1.0.0.
SkeelBerzins.jl is a universal solver for systems of onedimensional elliptic and parabolic nonlinear partial differential equations. The package provides an API similar to MATLAB’s pdepe and is based on the spatial discretization from the paper by Skeel and Berzins [1], which follows the method of lines together with a finiteelement discretization. The time discretization is done either by the DifferentialEquations.jl package or using the implicit Euler method (internal method).
The solver requires the PDE problem to be formulated as follows:
with associated boundary conditions expressed as:
for i=1, \; \cdots, \; npde (number of PDEs). Initial conditions must also be provided. For more information on the problem definition and assumptions, please consult the doc. For a list of various examples, check here.
Features of SkeelBerzins.jl:

the specific finiteelement discretization used ensures secondorder accuracy in space for various systems of coordinates (cartesian and polar). Convergence tests are available here.

wellintegrated into the SciML ecosystem with access to time integration methods through DifferentialEquations.jl and linear solvers via LinearSolve.jl.

designed with a focus on performance, using an allocationfree method for assembling the system of ODEs/DAEs.

interpolation available in time and space on the solution object (or by using the
pdeval
method). 
takes advantage of the operator’s sparsity to compute Jacobians in the Newton method and in the DifferentialEquations.jl solvers by using SparseDiffTools.jl.
Extended Capabilities:

possibility to solve pure elliptic problems by running one iteration of the implicit Euler method with a time step set to infinity.

variability of time discretizations through DifferentialEquations.jl.

Performance improvement, see figure below comparing SkeelBerzins.jl and MATLAB’s pdepe. Reproducible code is also available.

code designed to be compatible with general numeric datatypes.

possibility to use StaticArrays to define system of PDEs avoiding “heap allocations”.
Future extensions:

possiblity to define species only on a subregion of the general domain.

extension to solve twoscale problems.
Acknowledgments:
Many thanks to @jfu, for his valuable mentorship and numerous discussions.
Any suggestions for improving the package are more than welcome.
Best,
Grégoire Pourtier
Reference:
[1] Skeel, Robert D. and Berzins, Martin, “A Method for the Spatial Discretization of Parabolic Equations in One Space Variable”, SIAM Journal on Scientific and Statistical Computing, Vol. 11, 1990, pp.1–32.