I want to numerically integrate the equation $\partial_t u= a(t) \partial_xu+b\partial_{xxx}u+c$ to get $u(t)$. Is is preferable to use a difference formula of higher order of accuracy for spatial derivatives? Thanks. I am using a SemiDiscretization method called Method of Lines (vertical) to treat this problem.

You could use the method of lines to solve this PDE. If you use an explicit finite difference method, you will need to take a rather small time step (${\mathcal O(\Delta x^3))}$ due to the $u_{xxx}$ term. Given that you don't specify any boundary conditions, I will assume that you are solving the Cauchy problem. In that case (or in case of periodic boundary conditions), there is a much more efficient approach. Since your PDE is linear and the coefficients don't vary in space, you can solve in terms of Fourier modes. Decompose the initial data as usual: $$u(x,0) = \sum_k c_k\exp{ikx}.$$ Then use the ansatz $$u(x,0) = \sum_k g_k(t)\exp{ikx}$$ with $g_k(0) = c_k$. Substituting this in your PDE gives $$g_k'(t) = ik a(t)  ik^3b + c.$$ This can be solved very easily by numerical quadrature (or exactly, if $a(t)$ can be integrated symbolically). 

