Forward mode Jacobian

fwd.cpp

        for(int i = 0; i < y.size(); i++) {
          if(i > 0)
            y_var[i - 1].d_ = 0.0;
          y_var[i].d_ = 1.0;

          vector<fvar<double> > dy_dt_var = f_(t, y_var, thetavar, x_, x_int_, msgs_);
          if(i == 0) {
            for(int j = 0; j < y.size(); j++) {
              dy_dt(j) = dy_dt_var[j].val();
            }
          }

          for(int j = 0; j < y.size(); j++)
            Jy(j, i) = dy_dt_var[j].tangent();
        }

        y_var[y.size() - 1].d_ = 0.0;

        for(int i = 0; i < theta_.size(); i++) {
          if(i > 0)
            thetavar[i - 1].d_ = 0.0;
          thetavar[i].d_ = 1.0;

          vector<fvar<double> > dy = f_(t, y_var, thetavar, x_, x_int_, msgs_);

          for(int j = 0; j < dy.size(); j++)
            Jtheta(j, i) = dy[j].tangent();
        }

This goes in: stan/math/rev/mat/functor/ode_system.hpp in the:

template <typename Derived1, typename Derived2>
      inline void jacobian(double t, const std::vector<double>& y,
                           Eigen::MatrixBase<Derived1>& dy_dt,
                           Eigen::MatrixBase<Derived2>& Jy,
                           Eigen::MatrixBase<Derived2>& Jtheta)

function signature.