In this thesis we construct an accurate and effective numerical algorithm to solve the three classic 1D Hydraulic fracturing (HF) models: PKN, KGD and radial (also called penny-shaped). The solver works with power-law fluids (and Newtonian as a special case) and with impermeable rocks or with fluid leak-off modelled by classic Carter's formula. In KGD and radial models we also include the possibility to consider the effects of fluid induced shear stress on rock deformation and fracture propagation. The solver is adaptive in space and time, allowing at every stage of the computation to keep under control the error of the solution and to use the smallest possible number of approximation points. A brief description of the general HF problem is presented, and the three classic 1D HF models are described in detail. A simplified, time independent self-similar version of the 1D models is also considered. HF models are characterised by a moving boundary, for this reason we use a suitable normalisation of the variables to move the problem to a fixed interval. HF models are also characterised by an irregular (possibly singular), asymptotic behaviour of the solution at the crack mouth and at the crack tip, for this reason we use a carefully chosen change of variable to smooth the solution at the boundaries. Once the solution has been made smoother, it can be effectively approximated on the Chebyshev nodes with polynomials multiplied by the Jacobi weight function. Interpolation on Chebyshev nodes of a smooth function guarantees fast convergence and can be implemented efficiently using discrete cosine transform (or discrete Fourier transform). In the case of KGD and radial models we must also evaluate an integral operator that has an irregular kernel, this requires the use of another smoothing transformation. A fast and accurate way to evaluate the kernels of the integral operators using symmetric elliptic integrals is also proposed. The solution of the problem is further complicated by the fact that the asymptotic behaviour of the solution can change when passing from a regime to another: storage or leak-off dominated and toughness or viscosity dominated. Therefore, the function spaces of the approximants must be chosen even more carefully to be able to keep into account the behaviour of the solution at any regime. The self-similar and the time-dependent solvers that we propose are based on multigrid methods, where the solution is first found on a coarse grid and successively refined on denser grids, until the error is satisfactory. Time discretisation is done using implicit Runge Kutta methods, that allow high order of convergence while remaining stable in stiff problems. Time step strategy is also adaptive and step length is chosen dynamically depending on the solution. The solution is validated comparing it with some semi-analytical benchmarks present in literature and its convergence is thoroughly tested. Extensive computations prove that the numerical scheme is stable and efficient. It provides accurate results for all the three classic 1D hydraulic fracturing models, at all the regimes and with or without fluid leak-off. The algorithm can also be easily modified to work with different fluid and leak-off models. Finally, we run several simulations with oscillating pumping rate, oscillating leak-off coefficient and oscillating toughness, to see how this affects the behaviour of the solution compared to the case with constant parameters. In addition we investigate the effects of fluid induced shear stress on the fracture walls. We compare the results obtained with the modified KGD and radial models with those coming from the classic versions.