Variational methods for approximating disparity and optical flow are quite similar. Basically the same model that can be used for approximating optical flow with big displacements (late linearisation with warping) can also be used for approximating disparity. The only difference is that we don't expect vertical movement in the disparity case.
One complete chapter of my thesis is devoted to explaining how the variational optical flow, disparity and levelset models can be solved. In this chapter I show, step by step, how these models can be solved using different kinds of solvers (Jacobi, Gauss Seidel and Alternating Line Relaxation). Related to this chapter, I release the optical flow and disparity codes on this page...these codes are quite similar to the ones that I have used for generating the videos on this site. However, they are not 100% the same...there have been some slight modifications. For anyone interested in the actual code, I suggest first to have a look at the thesis in order to better understand the code itself! The codes are released under LGPL license. If you use the codes, please reference my papers...in my papers, on the other hand, I reference those papers upon which my work is based on. Thanx!!
The optical flow codes are as follows:
 Late linerisation opticalflow for large displacements.
 Early linearisation opticalflow method for small displacements (basically Horn&Schunck type of formulation).
 Late linearisation method for disparity.
Pseudo Code
Below I have included pseudo code of the latelinearization version of the optical flow, with warping, that can easily be modified for calculating stereo disparities.
 //
 //COARSETOFINE ALGORITHM FOR CALCULATING OPTICAL FLOW
 //Late linearisation (i.e. uses warping)
 //Robust error functions in both the data and the smoothness terms
 //Inputs are \(I_0\) and \(I_1\), number of scales \(scl\), and scaling factor \(sclFactor\)
 //
 \(\mathbf{INPUT:} \, I_0, \, I_1, \, scl, \, sclFactor \)
 \(\mathbf{OUTPUT:} \, (u,\, v) \)
 //Set \(u \) and \(v \) to zero
 \(u=0 \), \(v=0; \)
 //Create image pyramid
 \([ Iscl_0\{\} \, Iscl_1\{\} ] = pyramid(I_0, \, I_1, \, scl, \, sclFactor); \)
 //Coarsetofine loop
 WHILE( \(s=scl:1:1 \) )

 \(I_0 =Iscl_0\{s\} \), \(I_1 =Iscl_1\{s\}; \)
 //Warping loop
 WHILE( \( fstLoop \))

 //Warp image as per \(u \) and \(v \)
 \(I_{k,0}^w = warp( I_{k,0}, u, v ); \)
 Approximate derivatives for \(I_0^w \) and \(I_1 \): \(\dfrac{\partial I_k}{\partial t} = I_{k,1}  I_{k,0}^w \), \(\dfrac{\partial I_{k,0}^w}{\partial x} \), \(\dfrac{\partial I_{k,0}^w}{\partial y} \)
 //Reset \(du \) and \(dv \)
 \(du=0 \), \(dv=0 \)
 //Fixedpoint loop due to the robust error functions
 WHILE( \(sndLoop \))

 //Calculate penalizer function values for data
 \( \Psi{\prime} \Big( (E_k)_D \Big) \), where \(\left( E_k \right)_D = \left( \dfrac{ \partial I_k }{\partial t}  \dfrac{ \partial I_{k,0}^{w} }{\partial x}du  \dfrac{ \partial I_{k,0}^{w} }{\partial y}dv \right)^2 ; \)
 //Calculate the diffusion weights
 \( \Big[ \Psi{\prime} \big( E_R^{l,m} \big)_W \, \Psi{\prime} \big( E_R^{l,m} \big)_N \, \Psi{\prime} \big( E_R^{l,m} \big)_E \, \Psi{\prime} \big( E_R^{l,m} \big)_S \Big] = weights(u+du, v+dv); \)
 //Solve for new \(du \) and \(dv \)
 \( \begin{matrix} [du \, dv] = SOLVER ( & u, \, v, \, du, \, dv, \, nLoops, \, \dfrac{\partial I_k}{\partial t}, \, \dfrac{\partial I_{k,0}^w}{\partial x}, \, \dfrac{\partial I_{k,0}^w}{\partial y}, \, \Psi{\prime} \Big( (E_k)_D \Big), &\\ & \Psi{\prime} \big( E_R^{l,m} \big)_W, \, \Psi{\prime} \big( E_R^{l,m} \big)_N, &\\ & \Psi{\prime} \big( E_R^{l,m} \big)_S, \, \Psi{\prime} \big( E_R^{l,m} \big)_E & ); \end{matrix} \)
 ENDWHILE
 //Update \(u \) and \(v \)
 \(u=u+du \), \(v=v+dv; \)
 ENDWHILE
 //Interpolate (prolongate) solution
 IF( \(s1>0 \) )
 \([u \, v] = prolongate( u, \, v, \,sclFactor ); \)
 ENDIF
 ENDWHILE