This paper provides a theoretical foundation and numerical method for an inverse source problem of identifying a temporal load G(t) in a cantilever beam, with Neumann measured boundary output (i.e. bending moment) M(t) := -r(0)u(xx)(0, t), governed by the Euler-Bernoulli equation rho(x)u(tt)+ mu(x)u(t) + (r(x)u(xx) = F(x)G(t), (x, t) is an element of (0, l) x (0, T), subject to the initial u(x, 0) = u g (x, 0) = 0, x E (0,1), and the boundary conditions u(0, t) = u(x)(0, t) = 0, r(l)u(xx)(x,l) = (r(l)u(xx)(l,t)) x = 0, t is an element of [0, T] The spatial load F(x) is assumed to be known. The inverse problem is reformulated as an operator equation (Psi G)(t) = M(t), t is an element of (0, T), by introducing the input-output operator (Psi G)(t) := -r(0)u(xx)(0, t; G). 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 G(t) and the Neumann measured output M(t) 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 Frechet differentiability of the Tikhonov functional J(G) = (1 /2)parallel to r(0)u(xx)(0, .; G) M parallel to(2)(L2(0,T)). 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 Frechet 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.