Consider the fractional Brownian Motion (fBM) $B^H=\{B^H(t): t \in [0,1] \}$ with Hurst index $H\in (0,1)$. We construct a probability space supporting both $B^H$ and a fully simulatable process $\hat B_ε^H $ such that $$\sup_{t\in [0,1]}|B^H(t)-\hat B_ε^H(t)| \le ε$$ with probability one for any user specified error parameter $ε>0$. When $H>1/2$, we further enhance our error guarantee to the $α$-Hölder norm for any $α\in (1/2,H)$. This enables us to extend our algorithm to the simulation of fBM driven stochastic differential equations $Y=\{Y(t):t \in[0,1]\}$. Under mild regularity conditions on the drift and diffusion coefficients of $Y$, we construct a probability space supporting both $Y$ and a fully simulatable process $\hat Y_ε$ such that $$\sup_{t\in [0,1]}|Y(t)-\hat Y_ε^H(t)| \le ε$$ with probability one. Our algorithms enjoy the tolerance-enforcement feature, i.e., the error bounds can be updated sequentially. Thus, the algorithms can be readily combined with other simulation techniques like multilevel Monte Carlo to estimate expectation of functionals of fBMs efficiently.