@Carl Love 

Sorry for my careless. I have constructed a usable approximate initial solution.

The data given in initial solution is the values of variables, I pasted a picture bellow, it's a matrix form.

There are 70 columns because i used output=mesh in dsolve in a similar bvp problem . It is automaticly generate  70 different values of x. And I used its result as the initial approximate solution. 

Those x values one row of the file, and x-values include the two boundary points. 

My ODEs are all first order.

@Carl Love 

Thank you for your reply, but now I have some discrete data to be used for the  initial solution. For example,

there are 79 variables and the corresponding values are given.

I don't know how to construct a usable approximate initial solution. I have attached a short part of the data bellow.

Thank you very much for your reply!

initial_solution-5.txt

some_discrete_data_for_initial_solution.mw

@Rouben Rostamian  

thank you very much. I have seen this way to define a approximate initial solution. But my ODEs are very large, there are 73 unknown variables. So I think a discrete initial solution is more appropriate. Do you know how to define a discrete initial solution?

@Preben Alsholm thanks for your comment.

@Carl Love I tried to select several different parameters as the continuation parameter. But I failed to find one which works. And I observed if I change the value of one parameter I given, the code can run successfully. It's just a constant parameter, I have no idea why this situation would arise. Do you have any idea about this? 

@tomleslie 

sorry,I only saw the results you have pasted first time which display -0.1 for y1_15. I should check the code. 

Thank you very much for your help!

@acer 

thank you very much

@tomleslie 

it seems not right. In solnumeric2, lambda equal to 0.2, but in the solution y1_15 still equal to -0.1? it should equal to -0.2, doesn't it? I will try it tomorrow.

dsolve_code.mw

sorry for the code, I re-upload the original code.