diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C
index c4aad3c5..292d01a9 100644
--- a/src/SIM/SIMbase.C
+++ b/src/SIM/SIMbase.C
@@ -208,6 +208,9 @@ bool SIMbase::preprocess (const IntVec& ignored, bool fixDup)
   static int substep = 10;
   this->printHeading(substep);
 
+  if (mySol)
+    mySol->setupSecondarySolutions();
+
   // Perform some sub-class specific pre-preprocessing, if any
   this->preprocessA();
 
diff --git a/src/Utility/AnaSol.C b/src/Utility/AnaSol.C
index 558b3404..4f0cd71c 100644
--- a/src/Utility/AnaSol.C
+++ b/src/Utility/AnaSol.C
@@ -17,13 +17,51 @@
 #include "Utilities.h"
 #include "IFEM.h"
 #include "tinyxml.h"
+
+#include <autodiff/reverse/var.hpp>
+
 #ifdef HAS_HDF5
-#include "ProcessAdm.h"
-#include "HDF5Writer.h"
 #include "FieldFunctions.h"
 #endif
 
 
+namespace {
+
+//! \brief Function returning the gradient of a scalar or vector function.
+template<class FromFunc, class ToFunc, class Ret>
+class GradientFunc : public ToFunc
+{
+public:
+  //! \brief The constructor initializes the function to use.
+  GradientFunc(const FromFunc& f) : func(f) {}
+
+protected:
+  //! \brief Evaluates the gradient in a point.
+  Ret evaluate(const Vec3& X) const { return func.gradient(X); }
+
+private:
+  const FromFunc& func; //!< Reference to scalar/vector function
+};
+
+
+//! \brief Specialization for symmetric tensors.
+//! \details We need to construct a symmtensor from the returned tensor.
+template<>
+SymmTensor GradientFunc<VecFunc,STensorFunc,SymmTensor>::evaluate(const Vec3& X) const
+{
+  Tensor tmp = func.gradient(X);
+  const size_t nsd = tmp.dim();
+  SymmTensor ret(nsd);
+  for (size_t i = 1; i <= nsd; ++i)
+    for (size_t j = i; j <= nsd; ++j)
+      ret(i,j) = tmp(i,j);
+
+  return ret;
+}
+
+}
+
+
 AnaSol::AnaSol (RealFunc* s1, VecFunc* s2,
                 VecFunc* v1, TensorFunc* v2, STensorFunc* v3)
   : vecSol(v1), vecSecSol(v2), stressSol(v3)
@@ -92,8 +130,15 @@ AnaSol::AnaSol (const TiXmlElement* elem, bool scalarSol)
   const char* type = elem->Attribute("type");
   if (type && !strcasecmp(type,"fields"))
     this->parseFieldFunctions(elem,scalarSol);
-  else
-    this->parseExpressionFunctions(elem,scalarSol);
+  else {
+    bool useAD = false;
+    utl::getAttribute(elem, "autodiff", useAD);
+    utl::getAttribute(elem, "symmetric", symmetric);
+    if (useAD)
+      this->parseExpressionFunctions<autodiff::var>(elem,scalarSol);
+    else
+      this->parseExpressionFunctions<Real>(elem,scalarSol);
+  }
 }
 
 
@@ -128,8 +173,13 @@ void AnaSol::initPatch (size_t pIdx)
 }
 
 
+template<class Scalar>
 void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
 {
+  using EvalF = EvalFuncSpatial<Scalar>;
+  using VecF = EvalMultiFunction<VecFunc,Vec3,Scalar>;
+  using TensorF = EvalMultiFunction<TensorFunc,Tensor,Scalar>;
+  using STensorF = EvalMultiFunction<STensorFunc,SymmTensor,Scalar>;
   std::string variables;
   const TiXmlElement* var = elem->FirstChildElement("variables");
   if (var && var->FirstChild())
@@ -168,16 +218,16 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
     if (scalarSol)
     {
       if (type == "expression") {
-        scalSol.push_back(new EvalFunction((variables+primary).c_str()));
-        parseDerivatives(static_cast<EvalFunction*>(scalSol.back()),prim);
+        scalSol.push_back(new EvalF((variables+primary).c_str()));
+        parseDerivatives(static_cast<EvalF*>(scalSol.back()),prim);
       } else
         scalSol.push_back(utl::parseRealFunc(primary,type,false));
     }
     else
     {
       if (type == "expression") {
-        vecSol = new VecFuncExpr(primary,variables);
-        parseDerivatives(static_cast<VecFuncExpr*>(vecSol),prim);
+        vecSol = new VecF(primary,variables);
+        parseDerivatives(static_cast<VecF*>(vecSol),prim);
       } else
         vecSol = utl::parseVecFunc(primary, type);
     }
@@ -190,9 +240,9 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
     utl::getAttribute(prim, "type", type);
     std::string prType = (type == "expression" ? "" : "("+type+") ");
     std::string primary = prim->FirstChild()->Value();
-    IFEM::cout <<"\tScalar Primary " << prType << "=" << primary << std::endl;
+    IFEM::cout <<"\tScalar Primary" << prType << "=" << primary << std::endl;
     if (type == "expression")
-      scalSol.push_back(new EvalFunction((variables+primary).c_str()));
+      scalSol.push_back(new EvalF((variables+primary).c_str()));
     else
       scalSol.push_back(utl::parseRealFunc(primary, type, false));
     prim = prim->NextSiblingElement("scalarprimary");
@@ -209,18 +259,18 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
     if (scalarSol)
     {
       if (type == "expression") {
-        scalSecSol.push_back(new VecFuncExpr(secondary,variables));
-        parseDerivatives(static_cast<VecFuncExpr*>(scalSecSol.back()),sec);
+        scalSecSol.push_back(new VecF(secondary,variables));
+        parseDerivatives(static_cast<VecF*>(scalSecSol.back()),sec);
       } else
         scalSecSol.push_back(utl::parseVecFunc(secondary, type));
     }
     else
     {
-      if (type == "expression")
-        vecSecSol = new TensorFuncExpr(secondary,variables);
-      else
+      if (type == "expression") {
+        vecSecSol = new TensorF(secondary,variables);
+        parseDerivatives(static_cast<TensorF*>(vecSecSol),sec);
+      } else
         vecSecSol = utl::parseTensorFunc(secondary, type);
-      parseDerivatives(static_cast<TensorFuncExpr*>(vecSecSol),sec);
     }
   }
 
@@ -232,7 +282,7 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
     std::string secondary = sec->FirstChild()->Value();
     IFEM::cout <<"\tScalar Secondary="<< secondary << std::endl;
     if (type == "expression")
-      scalSecSol.push_back(new VecFuncExpr(secondary,variables));
+      scalSecSol.push_back(new VecF(secondary,variables));
     else
       scalSecSol.push_back(utl::parseVecFunc(secondary, type));
     sec = sec->NextSiblingElement("scalarsecondary");
@@ -241,10 +291,11 @@ void AnaSol::parseExpressionFunctions (const TiXmlElement* elem, bool scalarSol)
   const TiXmlElement* stress = elem->FirstChildElement("stress");
   if (stress && stress->FirstChild())
   {
+    symmetric = true;
     std::string sigma = stress->FirstChild()->Value();
     IFEM::cout <<"\tStress="<< sigma << std::endl;
-    stressSol = new STensorFuncExpr(sigma,variables);
-    parseDerivatives(static_cast<STensorFuncExpr*>(stressSol),stress);
+    stressSol = new STensorF(sigma,variables);
+    parseDerivatives(static_cast<STensorF*>(stressSol),stress);
   }
 }
 
@@ -310,3 +361,23 @@ void AnaSol::parseFieldFunctions (const TiXmlElement* elem, bool scalarSol)
             <<", no fields read."<< std::endl;
 #endif
 }
+
+
+void AnaSol::setupSecondarySolutions ()
+{
+  // if we are given a stress sol, no scalar secondaries should be registered.
+  // this is tailored to assumptions in elasticity applications.
+  if (!stressSol && !scalSol.empty()) {
+    scalSecSol.resize(scalSol.size());
+    for (size_t i = 0; i < scalSol.size(); ++i)
+      if (!scalSecSol[i])
+        scalSecSol[i] = new GradientFunc<RealFunc,VecFunc,Vec3>(*scalSol[i]);
+  }
+
+  if (vecSol) {
+    if (symmetric && !stressSol)
+      stressSol = new GradientFunc<VecFunc,STensorFunc,SymmTensor>(*vecSol);
+    else if (!symmetric && !vecSecSol)
+      vecSecSol = new GradientFunc<VecFunc,TensorFunc,Tensor>(*vecSol);
+  }
+}
diff --git a/src/Utility/AnaSol.h b/src/Utility/AnaSol.h
index ba5c1a10..b2489470 100644
--- a/src/Utility/AnaSol.h
+++ b/src/Utility/AnaSol.h
@@ -123,14 +123,22 @@ public:
   //! \brief Sets the patch to use.
   void initPatch(size_t pIdx);
 
+  //! \brief Make sure we have a secondary solution.
+  //! \details If none is given, we use derivation (automatic or finite difference)
+  //!          to obtain one.
+  virtual void setupSecondarySolutions();
+
 private:
   //! \brief Parses expression functions from XML definition.
+  template<class Scalar>
   void parseExpressionFunctions(const TiXmlElement* elem, bool scalarSol);
 
   //! \brief Parses field functions from XML definition.
   void parseFieldFunctions(const TiXmlElement* elem, bool scalarSol);
 
 protected:
+  bool symmetric = false; //!< True to use symmetric secondary solution
+
   std::vector<RealFunc*>   scalSol; //!< Primary scalar solution fields
   std::vector<VecFunc*> scalSecSol; //!< Secondary scalar solution fields
 
diff --git a/src/Utility/ExprFunctions.C b/src/Utility/ExprFunctions.C
index 3d3c91b9..990c2882 100644
--- a/src/Utility/ExprFunctions.C
+++ b/src/Utility/ExprFunctions.C
@@ -22,7 +22,7 @@
 #include <autodiff/reverse/var.hpp>
 
 
-template<> int EvalFunc::numError = 0;
+template<> int EvalFuncScalar<Real>::numError = 0; //!< Explicit instantiation
 
 
 namespace {
@@ -90,7 +90,7 @@ void ExprException (const ExprEval::Exception& exc, const char* task,
     std::cerr <<": Unknown exception";
   }
   std::cerr << std::endl;
-  EvalFunc::numError++;
+  EvalFuncScalar<Real>::numError++;
 }
 
 
diff --git a/src/Utility/ExprFunctions.h b/src/Utility/ExprFunctions.h
index 857b720e..59482cf9 100644
--- a/src/Utility/ExprFunctions.h
+++ b/src/Utility/ExprFunctions.h
@@ -99,14 +99,16 @@ class EvalFuncSpatial : public RealFunc
     Scalar* z; //!< Z-coordinate
     Scalar* t; //!< Time
 
+    //! \brief Returns a const ref to a member.
+    //! \param dir One-based index to member
     const Scalar& get(int dir) const
     {
       switch (dir) {
-        default:
-        case 1: return *x;
-        case 2: return *y;
-        case 3: return *z;
-        case 4: return *t;
+        case  1: return *x;
+        case  2: return *y;
+        case  3: return *z;
+        case  4: return *t;
+        default: return *x;
       }
     }
   };
@@ -268,4 +270,7 @@ using TensorFuncExpr = EvalMultiFunction<TensorFunc,Tensor,Real>;
 //! Symmetric tensor-valued function expression
 using STensorFuncExpr = EvalMultiFunction<STensorFunc,SymmTensor,Real>;
 
+//! \brief Explicit instantiation of error flag.
+template<> int EvalFunc::numError;
+
 #endif
diff --git a/src/Utility/Functions.C b/src/Utility/Functions.C
index 88cb5c71..dea2c4c4 100644
--- a/src/Utility/Functions.C
+++ b/src/Utility/Functions.C
@@ -717,7 +717,7 @@ static const ScalarFunc* parseFunction (const char* type, char* cline, Real C)
   {
     cline = strtok(cline,":");
     IFEM::cout << cline;
-    EvalFunc::numError = 0;
+    EvalFuncScalar<Real>::numError = 0;
     ScalarFunc* sf = new EvalFunc(cline,"t",C);
     if (EvalFunc::numError > 0)
     {
diff --git a/src/Utility/Test/TestAnaSol.C b/src/Utility/Test/TestAnaSol.C
index c9d8583c..92be542d 100644
--- a/src/Utility/Test/TestAnaSol.C
+++ b/src/Utility/Test/TestAnaSol.C
@@ -12,6 +12,8 @@
 
 #include "AnaSol.h"
 #include "Function.h"
+#include "Tensor.h"
+#include "TensorFunction.h"
 #include "Vec3Oper.h"
 #include "tinyxml.h"
 
@@ -48,15 +50,15 @@ TEST(TestAnaSol, ParseDerivatives)
   ASSERT_TRUE(v != nullptr);
 
   Vec3 X(1.0,2.0,3.0);
-  EXPECT_FLOAT_EQ((*v)(X),8.0);
-  EXPECT_FLOAT_EQ(v->deriv(X,1),3.0);
-  EXPECT_FLOAT_EQ(v->deriv(X,2),4.0);
-  EXPECT_FLOAT_EQ(v->deriv(X,3),1.0);
-  EXPECT_FLOAT_EQ(v->dderiv(X,1,1),6.0);
-  EXPECT_FLOAT_EQ(v->dderiv(X,2,2),2.0);
-  EXPECT_FLOAT_EQ(v->dderiv(X,3,3),0.0);
-  EXPECT_FLOAT_EQ(v->dderiv(X,1,2),0.0);
-  EXPECT_FLOAT_EQ(v->dderiv(X,3,2),0.0);
+  EXPECT_DOUBLE_EQ((*v)(X),8.0);
+  EXPECT_DOUBLE_EQ(v->deriv(X,1),3.0);
+  EXPECT_DOUBLE_EQ(v->deriv(X,2),4.0);
+  EXPECT_DOUBLE_EQ(v->deriv(X,3),1.0);
+  EXPECT_DOUBLE_EQ(v->dderiv(X,1,1),6.0);
+  EXPECT_DOUBLE_EQ(v->dderiv(X,2,2),2.0);
+  EXPECT_DOUBLE_EQ(v->dderiv(X,3,3),0.0);
+  EXPECT_DOUBLE_EQ(v->dderiv(X,1,2),0.0);
+  EXPECT_DOUBLE_EQ(v->dderiv(X,3,2),0.0);
 
   VecFunc* v2 = mySol.getScalarSecSol();
   ASSERT_TRUE(v2 != nullptr);
@@ -65,3 +67,124 @@ TEST(TestAnaSol, ParseDerivatives)
   EXPECT_TRUE(v2->deriv(X,2) == Vec3(0.0,2.0,0.0));
   EXPECT_TRUE(v2->dderiv(X,1,1) == Vec3(6.0,0.0,0.0));
 }
+
+
+TEST(TestAnaSol, ParseAD)
+{
+  const char* input = "<anasol type=\"expression\" autodiff=\"true\">"
+    "  <primary>"
+    "    x^3+y^2+z"
+    "  </primary>"
+    "</anasol>";
+
+  TiXmlDocument doc;
+  doc.Parse(input,nullptr,TIXML_ENCODING_UTF8);
+  const TiXmlElement* tag = doc.RootElement();
+  ASSERT_TRUE(tag != nullptr);
+  ASSERT_EQ(strcmp(tag->Value(),"anasol"),0);
+
+  AnaSol mySol(tag,true);
+  RealFunc* v = mySol.getScalarSol();
+  ASSERT_TRUE(v != nullptr);
+
+  mySol.setupSecondarySolutions();
+
+  Vec3 X(1.0,2.0,3.0);
+  const Vec3 grad1 = v->gradient(X);
+  const SymmTensor hess = v->hessian(X);
+  VecFunc* v2 = mySol.getScalarSecSol();
+  ASSERT_TRUE(v2 != nullptr);
+  const Vec3 grad2 = (*v2)(X);
+  EXPECT_DOUBLE_EQ((*v)(X),8.0);
+  EXPECT_DOUBLE_EQ(v->deriv(X,1),3.0);
+  EXPECT_DOUBLE_EQ(grad1[0],3.0);
+  EXPECT_DOUBLE_EQ(grad2[0],3.0);
+  EXPECT_DOUBLE_EQ(v->deriv(X,2),4.0);
+  EXPECT_DOUBLE_EQ(grad1[1],4.0);
+  EXPECT_DOUBLE_EQ(grad2[1],4.0);
+  EXPECT_DOUBLE_EQ(v->deriv(X,3),1.0);
+  EXPECT_DOUBLE_EQ(grad1[2],1.0);
+  EXPECT_DOUBLE_EQ(grad2[2],1.0);
+
+  EXPECT_DOUBLE_EQ(v->dderiv(X,1,1),6.0);
+  EXPECT_DOUBLE_EQ(hess(1,1), 6.0);
+  EXPECT_DOUBLE_EQ(v->dderiv(X,2,2),2.0);
+  EXPECT_DOUBLE_EQ(hess(2,2), 2.0);
+  EXPECT_DOUBLE_EQ(v->dderiv(X,3,3),0.0);
+  EXPECT_DOUBLE_EQ(hess(3,3), 0.0);
+  EXPECT_DOUBLE_EQ(v->dderiv(X,1,2),0.0);
+  EXPECT_DOUBLE_EQ(hess(1,2), 0.0);
+  EXPECT_DOUBLE_EQ(v->dderiv(X,3,2),0.0);
+  EXPECT_DOUBLE_EQ(hess(3,2), 0.0);
+}
+
+
+TEST(TestAnaSol, ParseFD)
+{
+  const char* input = "<anasol type=\"expression\">"
+    "  <primary>"
+    "    x^3+y^2+z"
+    "  </primary>"
+    "</anasol>";
+
+  TiXmlDocument doc;
+  doc.Parse(input,nullptr,TIXML_ENCODING_UTF8);
+  const TiXmlElement* tag = doc.RootElement();
+  ASSERT_TRUE(tag != nullptr);
+  ASSERT_EQ(strcmp(tag->Value(),"anasol"),0);
+
+  AnaSol mySol(tag,true);
+  RealFunc* v = mySol.getScalarSol();
+  ASSERT_TRUE(v != nullptr);
+
+  mySol.setupSecondarySolutions();
+
+  Vec3 X(1.0,2.0,3.0);
+  const Vec3 grad1 = v->gradient(X);
+  VecFunc* v2 = mySol.getScalarSecSol();
+  ASSERT_TRUE(v2 != nullptr);
+  const Vec3 grad2 = (*v2)(X);
+  EXPECT_DOUBLE_EQ((*v)(X),8.0);
+  EXPECT_NEAR(v->deriv(X,1),3.0,1e-6);
+  EXPECT_NEAR(grad1[0],3.0,1e-6);
+  EXPECT_NEAR(grad2[0],3.0,1e-6);
+  EXPECT_NEAR(v->deriv(X,2),4.0,1e-6);
+  EXPECT_NEAR(grad1[1],4.0,1e-6);
+  EXPECT_NEAR(grad2[1],4.0,1e-6);
+  EXPECT_NEAR(v->deriv(X,3),1.0,1e-6);
+  EXPECT_NEAR(grad1[2],1.0,1e-6);
+  EXPECT_NEAR(grad2[2],1.0,1e-6);
+}
+
+
+TEST(TestAnaSol, ParseADStress)
+{
+  const char* input = "<anasol type=\"expression\" symmetric=\"true\" autodiff=\"true\">"
+    "  <primary>"
+    "    x^3+y^2+z | x^2+y^3+z"
+    "  </primary>"
+    "</anasol>";
+
+  TiXmlDocument doc;
+  doc.Parse(input,nullptr,TIXML_ENCODING_UTF8);
+  const TiXmlElement* tag = doc.RootElement();
+  ASSERT_TRUE(tag != nullptr);
+  ASSERT_EQ(strcmp(tag->Value(),"anasol"),0);
+
+  AnaSol mySol(tag,false);
+  ASSERT_TRUE(mySol.getScalarSol() == nullptr);
+  ASSERT_TRUE(mySol.getVectorSol() != nullptr);
+
+  mySol.setupSecondarySolutions();
+
+  ASSERT_TRUE(mySol.getVectorSecSol() == nullptr);
+  const STensorFunc* v = mySol.getStressSol();
+  ASSERT_TRUE(v != nullptr);
+
+  Vec3 X(1.0,2.0,3.0);
+  const SymmTensor grad1 = (*v)(X);
+  EXPECT_DOUBLE_EQ(grad1(1,1), 3*X.x*X.x);
+  EXPECT_DOUBLE_EQ(grad1(1,2), 2*X.y);
+  EXPECT_DOUBLE_EQ(grad1(2,1), 2*X.y);
+  EXPECT_DOUBLE_EQ(grad1(2,2), 3*X.y*X.y);
+}