When developing efficient numerical methods for solving parabolic types of equations, severe temporal stability constraints on the time step are often required due to the high-order spatial derivatives and/or stiff reactions. The implicit integration factor (IIF) method, which treats spatial derivative terms explicitly and reaction terms implicitly, can provide excellent stability properties in time with nice accuracy. One major challenge for the IIF is the storage and calculation of the dense exponentials of the sparse discretization matrices restulted from the linear differential operators. The Krylov subspace approximations can overcome this shortcoming and greatly save computational cost and storage as long as the Arnoldi errors are smaller than the schemes local truncation errors. Most explicit numerical methods for solving fourth order partial differential equations (PDEs) have time constraints that require small time steps, which is not feasible for high dimensional problems. When dealing with time-dependent fourth order PDEs with linear high order terms and stiff lower order nonlinear terms, implicit methods are often used to solve these types of equations. Previously, implicit integration factor (IIF) methods are often a good choice as a class of efficient "exactly linear part" time discretization method for second order equations, namely advection-diffusion-reaction (ADR) equations. I provide linear analysis on the fully discretized fourth order equation using the implicit integration factor method and show the error analysis. Solving time dependent fourth order partial differential equations usually suffer from slow time evolution. Previous methods that solve fourth order equations explicitly need small time steps. I provide a solution to solving fourth order partial differential equations that have time steps of the same magnitude as the spatial grid size. . I provide an introduction to the various fields in which fourth order partial differential equations are used. I will give details of the implicit integration factor method that integrates fully nonlinear advection-diffusion-reaction equations exactly in time. I will present numerical results by solving 1-D Cahn-Hilliard equation and 1-D Kuramoto-Sivashinsky equation using both the second and third order implicit integration factor method (IIF). I also present several fourth order equations in 2-D and 3-D that have stiff reactions. I will discuss computational efficiency with Krylov subspace approximation and the benefits of having large time step sizes in all problems.