Thank you very much @mmcdara 5655. The boundary conditions have been included for u[1](y),w[1](y),theta[1](y), and C[1](y)

New.mw

Thanks @tomleslie. I actually make a type in the Energy equation which resulted in the error code. In anyways, thanks for your suggestion.

@tomleslie how can I overcome this, since it is inside the loop from the beginning. I'd appreciate if you put me through. Thanks

Thanks @tomleslie, it is now working.

Thank you @tomleslie 5388,

Maybe the problem is specifically with Maple 2018. because I can see here that its working perfectly but still not substituting in my Maple 2018. Any help please. 

@Carl Love Thank you very much for your analysis. As for the first error which leads to erroneous results is taken. But the other L[1] and L[2] were actually functions and not multiplications, The functions will operate on the equations to solve them but it's given out (table( [( 0 ) = proc (eta) .... which I can't decode. 

Please, help me check the to see if there is a way out.

Appreciate !!! 

@Carl Love Thanks. At this point what will you have be do because there is solutions for S[0]. I want it to substitute solutions of S[0] into S[1] and so forth. 

Can you plese help me out

Yes, According to the journal, psi was varied from -0.8 to 0.8 and x range from -0.2 to 0.7. While y = h.

Pls. any clue

@tomleslie 

I have paste it on your wall. Please have a look. Thankz

restart;
delta := proc (k) options operator, arrow; `if`(k = 0, 1, 0) end proc; Theta[0] := a; Theta[1] := 0; m := 15; beta := .1; S[h] := .1; N[c] := .1; P[e] := .1;
for k from 0 to 30 do Theta[k+2] := solve(add((delta(l)+beta*Theta[l])*(k-l+1)*(k-l+2)*Theta[k-l+2], l = 0 .. k)+beta*add((l+1)*Theta[l+1]*(k-l+1)*Theta[k-l+1], l = 0 .. k)-S[h]*add(Theta[l]*Theta[k-l], l = 0 .. k)-N[c]*Theta[k]-P[e]*(k+1)*Theta[k+1], Theta[k+2]) end do; theta := sum(Theta[j]*x^j, j = 0 .. m);
solve({eval(sum(Theta[j]*x^j, j = 0 .. m), x = 1) = 1}, [a]);
   [[a = 0.9156400273], [a = -7.939533630 + 6.634218100 I], 

     [a = -11.91470150 + 9.262617543 I], 

     [a = -6.095870030 + 3.866552541 I], 

     [a = -5.930605034 + 2.326548728 I], 

     [a = -21.61742707 + 7.420431340 I], 

     [a = -5.829137532 + 1.104465678 I], [a = -5.797547421], 

     [a = -13.45139065], [a = -5.829137532 - 1.104465678 I], 

     [a = -21.61742707 - 7.420431340 I], 

     [a = -5.930605034 - 2.326548728 I], 

     [a = -6.095870030 - 3.866552541 I], 

     [a = -11.91470150 - 9.262617543 I], 

     [a = -7.939533630 - 6.634218100 I]]
V1 := subs([a = .9156400273], theta);
                                    2                   3
      0.9156400273 + 0.08034511410 x  + 0.002453516662 x 

                           4                     5
         + 0.001497148773 x  + 0.00004119147538 x 

                             6                 -7  7
         + 0.00002199042947 x  + 7.624384961 10   x 

                         -7  8                 -9  9
         + 2.402068249 10   x  + 6.756951078 10   x 

                         -9  10                 -11  11
         + 1.805164817 10   x   + 3.430991798 10    x  

                         -12  12                 -13  13
         + 7.291252365 10    x   - 1.846619533 10    x  

                         -14  14                 -15  15
         - 3.605275446 10    x   - 4.279155395 10    x  
NULL;
restart;
delta := proc (k) options operator, arrow; `if`(k = 0, 1, 0) end proc; Theta[0] := a; Theta[1] := 0; m := 15; beta := .1; S[h] := .1; N[c] := .3; P[e] := .1;
for k from 0 to 30 do Theta[k+2] := solve(add((delta(l)+beta*Theta[l])*(k-l+1)*(k-l+2)*Theta[k-l+2], l = 0 .. k)+beta*add((l+1)*Theta[l+1]*(k-l+1)*Theta[k-l+1], l = 0 .. k)-S[h]*add(Theta[l]*Theta[k-l], l = 0 .. k)-N[c]*Theta[k]-P[e]*(k+1)*Theta[k+1], Theta[k+2]) end do; theta := sum(Theta[j]*x^j, j = 0 .. m);
solve({eval(sum(Theta[j]*x^j, j = 0 .. m), x = 1) = 1}, [a]);
  [[a = 0.8416571890], [a = -8.111465529 + 5.807801432 I], 

    [a = -11.58792285 + 8.168731482 I], 

    [a = -6.576261118 + 3.356629186 I], 

    [a = -20.14727468 + 6.567943736 I], 

    [a = -6.460202793 + 1.999658489 I], 

    [a = -6.391646400 + 0.9437066835 I], [a = -6.371259277], 

    [a = -13.50350873], [a = -6.391646400 - 0.9437066835 I], 

    [a = -6.460202793 - 1.999658489 I], 

    [a = -20.14727468 - 6.567943736 I], 

    [a = -6.576261118 - 3.356629186 I], 

    [a = -11.58792285 - 8.168731482 I], 

    [a = -8.111465529 - 5.807801432 I]]
V2 := subs([a = .8416571890], theta);
                             2                   3
0.8416571890 + 0.1491173506 x  + 0.004584703488 x 

                     4                    5                     6
   + 0.004448133301 x  + 0.0001180214600 x  + 0.00007207997076 x 

                        7                 -7  8
   + 0.000001662123933 x  + 8.410023876 10   x 

                   -8  9                 -9  10
   + 1.416870772 10   x  + 5.099405094 10   x  

                   -11  11                 -11  12
   - 6.845827501 10    x   - 1.899739485 10    x  

                   -12  13                 -13  14
   - 1.697500318 10    x   - 2.709937070 10    x  

                   -14  15
   + 3.503326091 10    x  
NULL;
restart;
delta := proc (k) options operator, arrow; `if`(k = 0, 1, 0) end proc; Theta[0] := a; Theta[1] := 0; m := 15; beta := .1; S[h] := .1; N[c] := 1; P[e] := .1;
for k from 0 to 30 do Theta[k+2] := solve(add((delta(l)+beta*Theta[l])*(k-l+1)*(k-l+2)*Theta[k-l+2], l = 0 .. k)+beta*add((l+1)*Theta[l+1]*(k-l+1)*Theta[k-l+1], l = 0 .. k)-S[h]*add(Theta[l]*Theta[k-l], l = 0 .. k)-N[c]*Theta[k]-P[e]*(k+1)*Theta[k+1], Theta[k+2]) end do; theta := sum(Theta[j]*x^j, j = 0 .. m);
solve({eval(sum(Theta[j]*x^j, j = 0 .. m), x = 1) = 1}, [a]);
  [[a = 0.6433384619], [a = -8.510805022 + 3.728372760 I], 

    [a = -12.93300935 + 3.848351148 I], 

    [a = -6.762161798 + 1.488424796 I], 

    [a = -10.02203491 + 0.2715688502 I], 

    [a = -10.27720548 + 0.1400907567 I], 

    [a = -9.790842845 + 0.1299464236 I], [a = -13.71571191], 

    [a = -9.790842845 - 0.1299464236 I], 

    [a = -10.27720548 - 0.1400907567 I], 

    [a = -10.02203491 - 0.2715688502 I], 

    [a = -6.762161798 - 1.488424796 I], 

    [a = -12.93300935 - 3.848351148 I], 

    [a = -8.510805022 - 3.728372760 I]]
V3 := subs([a = .6433384619], theta);
                             2                  3
0.6433384619 + 0.3216692310 x  + 0.01007419593 x 

                    4                    5                    6
   + 0.02380184188 x  + 0.0006769510602 x  + 0.0004518870434 x 

                        7                      8
   - 0.000005332330250 x  - 0.000006556130624 x 

                   -7  9                 -9  10
   - 7.333015310 10   x  - 3.879060316 10   x  

                   -8  11                 -8  12
   + 3.682027989 10   x   + 1.903529658 10   x  

                   -9  13                 -10  14
   + 1.413925023 10   x   - 2.257537367 10    x  

                   -10  15
   - 1.204519535 10    x  
NULL;
restart;
delta := proc (k) options operator, arrow; `if`(k = 0, 1, 0) end proc; Theta[0] := a; Theta[1] := 0; m := 15; beta := .1; S[h] := .1; N[c] := 5; P[e] := .1;
for k from 0 to 30 do Theta[k+2] := solve(add((delta(l)+beta*Theta[l])*(k-l+1)*(k-l+2)*Theta[k-l+2], l = 0 .. k)+beta*add((l+1)*Theta[l+1]*(k-l+1)*Theta[k-l+1], l = 0 .. k)-S[h]*add(Theta[l]*Theta[k-l], l = 0 .. k)-N[c]*Theta[k]-P[e]*(k+1)*Theta[k+1], Theta[k+2]) end do; theta := sum(Theta[j]*x^j, j = 0 .. m);
solve({eval(sum(Theta[j]*x^j, j = 0 .. m), x = 1) = 1}, [a]);
   [[a = 0.2160419278], [a = 0.9356057884 + 7.937063751 I], 

     [a = -4.929433007 + 20.39139370 I], 

     [a = -0.8369403204 + 3.385911692 I], 

     [a = -20.56438138 + 17.40656490 I], 

     [a = -1.554038706 + 1.218109435 I], 

     [a = -26.85780794 + 8.520367031 I], [a = -13.15066763], 

     [a = -31.01207367], [a = -26.85780794 - 8.520367031 I], 

     [a = -1.554038706 - 1.218109435 I], 

     [a = -20.56438138 - 17.40656490 I], 

     [a = -0.8369403204 - 3.385911692 I], 

     [a = -4.929433007 - 20.39139370 I], 

     [a = 0.9356057884 - 7.937063751 I]]
V4 := subs([a = .2160419278], theta);
                             2                  3
0.2160419278 + 0.5309674030 x  + 0.01732462882 x 

                   4                   5                  6
   + 0.2050550199 x  + 0.007390146595 x  + 0.02411044982 x 

                      7                    8                    9
   + 0.0005168127670 x  - 0.0008106869795 x  - 0.0001691360028 x 

                      10                     11
   - 0.0004176218493 x   - 0.00002502709502 x  

                   -7  12                      13
   + 3.647584125 10   x   + 0.000004207621358 x  

                       14                      15
   + 0.00001036406296 x   + 0.000001146428419 x  
NULL;
restart;
delta := proc (k) options operator, arrow; `if`(k = 0, 1, 0) end proc; Theta[0] := a; Theta[1] := 0; m := 15; beta := .1; S[h] := .1; N[c] := 10; P[e] := .1;
for k from 0 to 30 do Theta[k+2] := solve(add((delta(l)+beta*Theta[l])*(k-l+1)*(k-l+2)*Theta[k-l+2], l = 0 .. k)+beta*add((l+1)*Theta[l+1]*(k-l+1)*Theta[k-l+1], l = 0 .. k)-S[h]*add(Theta[l]*Theta[k-l], l = 0 .. k)-N[c]*Theta[k]-P[e]*(k+1)*Theta[k+1], Theta[k+2]) end do; theta := sum(Theta[j]*x^j, j = 0 .. m);
solve({eval(sum(Theta[j]*x^j, j = 0 .. m), x = 1) = 1}, [a]);
   [[a = 0.08667607473], [a = 4.512623776 + 3.614155013 I], 

     [a = 15.19532498 + 35.78098906 I], 

     [a = 0.7041394050 + 2.079680542 I], 

     [a = -0.5265761857 + 0.8179641462 I], 

     [a = -26.17482452 + 37.96400497 I], 

     [a = -43.16166445 + 18.91197981 I], [a = -13.18794443], 

     [a = -54.42771599], [a = -43.16166445 - 18.91197981 I], 

     [a = -26.17482452 - 37.96400497 I], 

     [a = -0.5265761857 - 0.8179641462 I], 

     [a = 0.7041394050 - 2.079680542 I], 

     [a = 15.19532498 - 35.78098906 I], 

     [a = 4.512623776 - 3.614155013 I]]
V5 := subs([a = 0.8667607473e-1], theta);
                              2                  3
0.08667607473 + 0.4300286912 x  + 0.01421111337 x 

                   4                  5                 6
   + 0.3470791612 x  + 0.01333277563 x  + 0.1009219890 x 

                     7                   8                    9
   + 0.003553431508 x  + 0.008177891042 x  - 0.0001428631684 x 

                     10                    11
   - 0.002705439006 x   - 0.0002847049854 x  

                      12                     13
   - 0.0008232307331 x   - 0.00004655348071 x  

                        14                     15
   + 0.000008501910885 x   + 0.00001316315336 x  
NULL;

plot({V1, V2, V3, V4, V5}, x = 0 .. 1, style = point, symbol = "asterisk", colour = ["blue", "red", "black", "purple", "green"]);
Warning, expecting only range variable x in expression V1 to be plotted but found name V1

@Kitonum 

Yes, I know it is a series of cos(parameters). But, I need to convert to the form "sine or cosine" if you have any idea.

Thanx

@Kitonum, Thank you for your suggestions. But the problem is the differential transform of Laguerre's equation in which the delta(i-1) can not be replace by delta(i). I have computed it manually but i need to know how to code it in maple. Possibly there is anoter way out please.

@tomleslie. Thank you very much for your suggestions. That was exactly what happened 

Thank U very much. It really worked

@nidojan, I'm very grateful. I think I will find a clue from this.