I just upgraded to Mathematica 13.1 and encountered the following issue (this is just code from the first few lines of the FEM programming documentation):


{state} = 
 NDSolve`ProcessEquations[{Laplacian[u[x, y], {x, y}] == 1, 
   DirichletCondition[u[x, y] == 0, True]}, u, {x, 0, 1}, {y, 0, 1}]
(* {NDSolve`StateData["<" "SteadyState" ">"]} *)

femdata = state["FiniteElementData"]
(* FiniteElementData["<" 1281 ">"] *)


state["FiniteElementData"]["Solution"] // Short

which resulted in the following error

During evaluation of In[16]:= FiniteElementData::nomthd: There is no method Solution for FiniteElementData objects.

Where did the solution data go? Is this a known issue?

  • 4
    $\begingroup$ The "Solution" method is already missing at least since v12.3. Not sure if it's a bug or compatibility issue. If I guess it right, the solution data can still be extracted with state["SolutionData"][[1, 3]] or state["SolutionVector"][[1]]. $\endgroup$
    – xzczd
    Commented Oct 26, 2022 at 7:53
  • $\begingroup$ I see! I did not see those listed when I executed state["Properties"]. I guess these are hidden. Thanks! $\endgroup$
    – Will.Mo
    Commented Oct 26, 2022 at 8:08

1 Answer 1


In essence this is a documentation issue. To extract the solution, there are two scenarios. Either you have a NDSolve state object (like you have shown) or you use a lower level function like PDESolve.

Consider this case:

{state} = 
  NDSolve`ProcessEquations[{Laplacian[u[x, y], {x, y}] == 1, 
    DirichletCondition[u[x, y] == 0, True]}, 
   u, {x, 0, 1}, {y, 0, 1}];

On this level you should use the solution data from the state object.

{forwardSolution} = state["SolutionData"];

For more information on the solution data see the documentation of the solution data object.

The actual data extraction should then happen with SolutionDataComponent

NDSolve`SolutionDataComponent[forwardSolution, "DependentVariables"]

The second case is when you use low level functions, like in this case:

nr = ToNumericalRegion[Rectangle[{0, 0}, {1, 1/2}]];
vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData["Space" -> nr];
methodData = InitializePDEMethodData[vd, sd]
pdeData = InitializePDECoefficients[vd, sd, 
   "DiffusionCoefficients" -> {{-IdentityMatrix[2]}}];
bcData = InitializeBoundaryConditions[vd, 
   sd, {{DirichletCondition[u[x, y] == 0., x == 0], 
     DirichletCondition[u[x, y] == 1., x == 1]}}];

Now we call PDESolve:

sdNew = PDESolve[pdeData, bcData, vd, sd, methodData];

Note that PDESolve returns a new solution data object.

sdNew === sd
(* False *)

You can then extract the solution vector from the object in the same manner as above:

NDSolve`SolutionDataComponent[sdNew, "DependentVariables"] // Short

This is a much cleaner process than what we had before. I have updated the reference pages to reflect this.


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.