Simulating a Hamilton-Jacobi-Bellman Say I have solved an HJB of the form:

$\rho V(k) = \max_c g(c) + V'(k)(z - c)$

I have calibrated $\rho$ to monthly parameters. I would like to simulate the development of $k$. I start with $k(0)$. However, unlike in discrete-time, I'm not sure what happens next.

Is $k(1) = k(0) + (z-c(0))$, with $k(1)$ being the value of $k$ for the next month? Or do I need to somehow integrate along the lines of $k(1) = \int_0^1 (z-c(t))dt + k(0)$? There's many references on solving these things, but I couldn't find any references on simulating it.
Expanding on jmbejara's answer and your question about integration, there are really three steps here:


Solve the HJB equation $\rho V(k) = \max_c g(c)+V'(k)(z-c)$. (You're already done this.)
Obtain the policy function $c(k)=\text{argmax}_c g(c)+V'(k)(z-c)$. (You presumably already did this while solving the HJB.)
Substitute $c(k)$ into the law of motion for $k$ and numerically solve that differential equation for $k(t)$:$$\dot{k}=z-c(k)$$


It sounds like you already had (1) and (2) and were just unsure about what exactly to do in (3). If you're working in a deterministic environment (and from the HJB it looks like you are), no discretization on your part is necessary: you can just use a canned ODE solver from Matlab or any other software package to numerically solve the differential equation above for $k(t)$.

To be a little more specific: 


If the HJB is simple and you have an analytical solution for $V(k)$, and thus can easily compute $c(k)$ as well, that's great. You can just plug $c(k)$ into the law of motion for $k$ and tell your software package of choice to solve it and obtain the trajectory of $k(t)$.
If you had to numerically solve the HJB, then that's great because it means that you probably have already used some kind of numerical ODE solution technique and will be familiar with it. On the minus side, it means that when you initially obtain $V(k)$ and $c(k)$, you'll only get it for a grid of $k$, rather than having some kind of formula you can apply everywhere - which is a problem, because you'll need to be able to determine $c(k)$ at any point $k$ in a continuum if you're going to numerically solve the $\dot{k}$ differential equation for $k(t)$. The simplest way around this, to my knowledge, is to fit a spline to (or use some other approximation method on) the discrete grid of $(k_{grid},c(k_{grid}))$ values that you have. Then you can obtain an approximate value of $c(k)$ anywhere that's covered by the spline, and it will be extremely accurate assuming that your ODE is well-behaved and you have a dense grid.


Things get more complicated in a stochastic environment, but your case looks deterministic. (And to be clear, I'm just a dabbler here: numerical methods are not my area of expertise. But this is what I know.)