Motivated from an empirical analysis of data collected by a smoking cession study,we propose a joint model (JM) of survival data and multiple longitudinal covariate pro-cesses, develop an estimation procedure for this model using likelihood-based approach,and further establish the consistency and asymptotic normality of the resulting estimate. Computation for the proposed likelihood-based approach in the joint modeling is particularly challenging since the estimation procedure involves numerical integration over multi-dimensional space for the random effects in the JM. Existing numerical inte-gration methods become ineffective or infeasible for the JM. We introduce a numerical integration method based on computer experimental designs for the JM. We conduct Monte Carlo simulation to examine the finite sample performance of the procedure and compare the new numerical integration method with existing ones. We further illustrate the proposed procedure via an empirical study of smoking cession data.