Source identification problems in a system governed by Euler-Bernoulli beam equation rho(x)u(tt) + (r(x)u(xx))(xx) = F(x)H(t), (x, t) is an element of (0, l) x (0, T), from available boundary observation (measured data), namely, from measured slope theta(t) := u(x)(0, t) at x = 0, are considered. We propose a new approach to identifying the unknown temporal (H(t)) and spatial (F(x)) loads. This novel approach is based on weak solution theory for PDEs and quasi-solution method for inverse problems combined with the adjoint method. It allows to construct not only a mathematical theory of inverse source problems for Euler-Bernoulli beam, but also an effective numerical algorithm for reconstruction of unknown loads. Introducing the input-output operators (Phi H) := u(x)(0, t; H) and (Psi F)(t) := u(x)(0, t; F), t is an element of (0, T), we show that both operators are compact. Based on this result and general regularization theory, we prove an existence of unique solutions of the regularized normal equations (Phi Phi* + alpha l)H-alpha = Phi*theta and (Psi*Psi + alpha l)F-alpha = Psi*theta. Then we develop the adjoint problem approach to prove Frechet differentiability of the corresponding cost functionals and Lipschitz continuity of the Frechet gradients. Derived explicit gradient formulas via the adjoint problem solution and known load, allow use of gradient type convergent iterative algorithms. Results of numerical simulations for benchmark problems illustrate robustness and high accuracy of the algorithm based on the proposed approach. (C) 2016 Elsevier Ltd. All rights reserved.