Could you please advise some literature about the numerical method of lines (MOL) for parabolic PDEs? It is a method of solving PDEs with discretizing only by space but not by time. A system of ODEs is obtained and it can be solved with, for example, Runge-Kutta method.
Do you know some solvers (programs, libraries) for applying the method of lines? Is this method applicable to multi-dimensional (2D, 3D) problems?
A nice reference for this is Chapter 10 of LeVeque's finite difference book. Of course, it only covers basic finite difference approaches, and there are plenty of others (all within the method of lines framework). The method of lines is indeed applicable in multiple dimensions, and two-dimensional problems are discussed in the reference just given.
Most of the methods and libraries that exist are based on the method of lines, so if you want meaningful recommendations you'll need to ask a more specific question.
This is a fairly big question. The basic idea is that in solving $\partial_t u = \mathcal{L}u$, you approximate $\mathcal{L}$ with a matrix $A$, and solve $$ u' = Au. $$ In this formulation you have two main issues:
These issues are not at all unique to the method of lines, or to parabolic equations. Books like "Solving ordinary differential equations" by Hairer and Wanner will cover this; "Finite difference methods" by LeVeque; there's an unfinished textbook available for free online here by Trefethen; other PDE textbooks cover this as well.
Note that just an explicit Runge-Kutta method to unlikely to work directly: the system of ODEs is typically stiff, sparse, and very large.
The question of how to solve the resulting systems of linear equations, roughly, can be addressed in different ways. If you can solve $A^{-1}x$ directly and efficiently, there is no problem. There are (diverse!) iterative methods, which work well if you can find a good preconditioner for the matrix. There are also different splitting methods, where you can write $A=A_1+A_2+\cdots$, and define a numerical solution to the ODE that depends only on solutions of systems $A_k^{-1}x$. You can pick $A_k$ to be convenient to solve, such as the individual parts of the differential operator going along individual dimensions. Alternating direction implicit methods are an example of this approach.
For software, most PDE packages do something like this, something like CVODE (in SUNDIALS) might be appropriate.
EDIT For nonlinear systems, the procedure is roughly the same, only instead of needing to solve systems like $A^{-1}x$ you need to solve systems of nonlinear equations with a Newton-like method, where the solution of $J^{-1}x$ with the Jacobian matrix needs to be efficient; again the Jacobian is likely to be sparse and very large. In this context, CVODE looks appropriate.
EDIT As David Ketcheson correctly points out, my answer assumes many things that are not explicitly specified in the question, which would make this answer quite misleading if those assumptions don't hold.
I don't have enough reputation to comment, but I hope this answer helps you. Personally, I find that I learn best through example. Try looking at the code here to see how MOL was implemented in Python with centered finite difference approximation (an ODE solver was used).
If you'd like to learn more about solving ODEs, parabolic PDEs, and MOL, I suggest Heath's Scientific Computation (Chapter 11) which also provides some examples on solving MOL. Good luck!
