title: Solve an intergral in R: the integrand is a function of the solution of an ODE
body: <p>I would like to compute an integral where the integrand is a function of the solution of an ODE. </p>
  
  <p>In order to solve the integral, R needs to solve an ODE for each value the integration algorithm uses. This is what I have done so far:</p>
  
  <pre><code>require(deSolve)
  
  # Function to be passed to zvode in order to solve the ODE  
  ODESR &lt;- function(t, state, parameters) {
     with(as.list(c(state, parameters)),{
         dPSI &lt;- -kappa*PSI+0.5*sigma^2*PSI^2
         dPHI &lt;- kappa*theta*PSI
         return(list(c(dPSI, dPHI)))
     })
  }
  
  # For a given value of p this code should return the solution of the integral
  pdfSRP &lt;- function (p) {
     integrand &lt;- function (u) {
         state &lt;- c(PSI = u*1i, PHI = 0)
         out &lt;- as.complex(zvode(y = state, times = times, parms = parameters, fun = ODESR)[2, 2:3])
         Re(exp(out[2] + out[1]*x)*exp(-u*1i*p))
     }
  integrate(f = integrand, lower = -Inf, upper = Inf)$value/(2*pi)
  }
  </code></pre>
  
  <p>For the following given values:</p>
  
  <pre><code>parameters &lt;- c(kappa = 1, theta = 0.035, sigma = 0.05)
  times &lt;- c(0,1)
  x &lt;- 0.1
  </code></pre>
  
  <p>running:</p>
  
  <pre><code>pdfSRP(p = 2)
  </code></pre>
  
  <p>produces the following error:</p>
  
  <pre><code>Error in eval(expr, envir, enclos) : object 'PSI' not found
  </code></pre>
  
  <p>I just cannot figure out why. I'm quite sure it is due to a syntax error, because running:</p>
  
  <pre><code>integrand &lt;- function (u) {
         state &lt;- c(PSI = u*1i, PHI = 0)
         out &lt;- as.complex(zvode(y = state, times = times, parms = parameters, fun = ODESR)[2, 2:3])
         Re(exp(out[2] + out[1]*x)*exp(-u*1i*p))
  }
  </code></pre>
  
  <p>with <code>p &lt;- 2</code> and (for example) <code>u &lt;- 3</code> works.</p>
  
  <p>Can you help me spot the mistake?</p>