I will present an algorithm allowing for a FULLY AUTOMATIC solution of elliptic and (time-harmonic) Maxwell’s equations with the hp-adaptive edge finite elements. With no interaction from the user, except for defining the problem, the algorithm delivers a sequence of optimally refined hp meshes, aimed at delivering a high quality solution (typically one percent error measured in energy norm, relative to the energy of the solution), with a minimum number of degrees-of-freedom. Construction of the meshes involves not only identifying which elements to refine, but also deciding whether to increase the order of approximation p, or break the elements, decreasing element size h. Additionally, one has to decide between ISOTROPIC or ANISOTROPIC refinements, an issue critical in capturing effectively various types of singularities and boundary layers. The algorithm (to the best of our knowledge - first of its kind), is based on identifying the optimal refinements, by minimizing the corresponding hp-interpolation error for a reference solution which is obtained by solving the problem on a mesh obtained by globally refining the current mesh in BOTH h and p. Critical to the efficiency of the method is the use of a two-grid solver for obtaining the reference solution. The method has been successfully tested on a series of 2D test problems, delivering (known from elsewhere) optimal meshes, and exponential convergence rates. The test problems to be presented include: a/ for H1-conforming elements: - 2D Laplace (Poisson) equation, - 2D non-homogeneous, highly anisotropic heat conduction, - 2D system of linear elaticity equations, - 3D axisymmetric problem for Maxwell equations, reduced to a 2D elliptic problem, b/ for H(curl)-conforming elements: - 2D diffraction problem for Maxwell equations. Additionally, I will show results on combining hp- adaptivity with a goal-oriented adaptive strategy, and (hopefully) first 3D results for 3D Fichera’s corner problem.