This paper provides a theoretical foundation and numerical method for an inverse source problem of identifying a temporal load in a cantilever beam, with Neumann measured boundary output (i.e. bending moment) , governed by the Euler–Bernoulli equation , , subject to the initial u(x, 0) = ut(x, 0) = 0, , and the boundary conditions u(0, t) = ux(0, t) = 0, , . The spatial load is assumed to be known. The inverse problem is reformulated as an operator equation , , by introducing the input–output operator . Our approach is based on detailed analysis of this operator by developing the regular weak solution theory for the direct problem. We show how to choose the temporal load and the Neumann measured output in order to guarantee solvability of the inverse problem. We introduce a weaker solution to the adjoint problem corresponding to the inverse problem and prove Fréchet differentiability of the Tikhonov functional . An explicit gradient formula for the Tikhonov functional is derived by making use of the unique weaker solution of the adjoint problem. Furthermore, under additional regularity and consistency conditions, the Lipschitz continuity of the Fréchet gradient is proved. This property allows use of the conjugate gradient algorithm (CGA) which is the main computational tool in solving inverse problems. A numerical method based on the finite-element discretization and the CGA is developed for the solution of the inverse source problem. Numerical examples with random noisy measured output are presented to illustrate the validity and effectiveness of the proposed approach.