Hello!
As a first step I'd estimate whether you actually need to re-use the matrix from the previous step. Maybe, you'll get results quickly enough when simply re-computing factorization at each step (= for each matrix) and there is no need to complicate things.
The question is then what remains the same in your problem when the matrix size changes. Do you expect previous matrix to be a good approximation of the current one? Then possibly you can use the factorization once you've computed it as a preconditioner for the next problem but you'll need a transition matrix which will match you meshes at different time steps (I imagine your problem size changes due to the mesh changing). In this case, you'll have your favorite iterative method (say, CG or BiCGStab) with a preconditioner which will involve calling PARDISO when computing the preconditioner action.
Also, if you matrix changes are due to the localized changes in the mesh, you can do the same, or exploit the Schur complement (re-compute only a smaller part of the factorization). See https://software.intel.com/en-us/articles/intel-mkl-support-to-new-functionality-schur-complement
Hope this is useful.
Best,
Kirill
