Thanks for your advice. I am now trying to use NDSolve. It's indeed much faster than using DSolve. But I've encountered another problem. For 6 second-order ODEs with 12 boundary conditions, I can obtain the exact symbolic solutions in relatively less time by using DSolve. Then I plotted the first derivative of the symbolic solutions, from which I got the right graph. But when I used NDSolve to solve the same problem, I got nothing after I plotted the first derivative of the numerical solutions. Why so? Besides, the graph of the symbolic solutions is the same as the graph of the numerical solutions. Then why the graphs of the first derivative of the two solutions are totally different? Looking forward to your reply. Many thanks!!