At present I am writing a controller by locally linearising the robot model using finite differences. \dot{x} = f(x,u) is locally linearised to \dot{x} = Ax + Bu. Is there a better way to get A and B or the dynamic model from mujoco (this is for optimization based approaches)?
Check out the new sections in the Programming Guide in the MuJoCo Pro chapter -- there is a code sample called derivative.cpp and detailed explanation on how to compute finite differences efficiently using multi-threading and controlled warm-starts of the iterative solver. Sometime this summer there will be a release with optimization-related features, including analytical computation of the B matrix and parts of the A matrix (some finite differences will still be needed at the end, because differentiating the collision detectors analytically is not feasible). In general, MuJoCo 2.x will be mostly about model-based optimization, although the underlying simulation will still continue to be improved.