@tomleslie also the points obtained from BCs and ICs are not in the matrices which you generate.

@tomleslie dear thanks for your positive response. for 1/(h=5) and 1/(k=5) the program give good results but for higher values (arbitrary values) of h and k results seems not good please see your sent program.

@tomleslie thanks for your answer. Now I want to minimize the size of 3D plot. I know how to do it for 2D plot, but don't know in 3D. Please help me to do this. I am waiting your response.

@Carl Love yes I also need legend please guide me how to do. Give nine curve arbitrary name.

Dear @Carl Love thanks for your deep answer about alpha. I got it and do the changing as per your suggestion. But you said you solve the system for N=100 in on mint. Can you share me that maple? I want to know who computational work reduces. 

Dear @Markiyan Hirnyk first of all I am very thankful your positive response. I did the changing according your instruction but could'nt find results becuase maple sheet take alot of time to evalute. Please see the attachment and see why they takes to much time. I want to produce reuslts for M=N=100. Again thanks in advance.

ImprovedEquations.mw

Dear @Carl Love first of all I am very thankful your positive response. I did the changing according your instruction but could'nt find results becuase maple sheet take alot of time to evalute. Please see the attachment and see why they takes to much time. I want to produce reuslts for M=N=100. Again thanks in advance.

Equations.mw

Dear @Kitonum first of all thanks alot for your positive answer. You almost fix my problem but I need answer in form that the denominator seprate with each term then simple.

dear @taro first of all thanks for your answer! I want to covert diff(u,y) into diff(U,Y) with u=v/d^2*U(Y), and y=d*Y. Manually I do in the attachment, which can help to fix the problem. 

manual.pdf

I forget to attach the file please see now
Code_paper_1.mw

Dear @Carl Love thanks for your answer. I need the inverse laplace of F(q, x, alpha) that is invlaplace(F(q, x, alpha),q,t) then I want to plot it for fix values of t and alpha for x=0..1. 

Dear @vv; thanks for your precious time and valuable suggestions. 

I know each and everything. If you can solve my problem then reply. I am not free to see your bloody answer.

@Preben Alsholm 

Thanks for your kind answer. Now my one more question is which method (buiton command) is best to solve system of nonlinear ODEs numerically. Beacuse some time I face the error like this mehtod is not converges.

@tomleslie 

Dear Sir! I want to see the values of F(eta), G(eta) and H(eta) form 0 to 1 with 0.1 difference. I need the accuracy till 13 decimal places please correct the program accordingly, Thanks for your psoitive response.