@jrive 

I would use Syrup's "ladder syntax" as a general rule. It is (ususally) simpler to produce a netlist just by specifying components in the form

componentName nodeName1 nodeName2 componentValue

for example

R1 1 2 10

The following code (and the attached syrup2.mw) is about as simple as I can make it (although it does show that producing plots *may* be simpler using the DynamicSystems:-BodePlot() command, which saves typing a lot of options into a basic plot() command

  with(Syrup):
  with(plots):
  ckt:="
         V1 1 0 1
         R1 1 3 50
         L1 3 2 .96e-5
         Cp 2 0 .83e-10
         L2 2 5 500e-6
         R2 5 4 0.2
         RL 4 6 1.343
         LL 6 0 0.155e-05
         .end
       ":
  sol:= Solve(ckt, ac):
#
# Extract the "output node from the above solution
# assuming it is the one numbered "4"
#
  g1:= eval( v[4], sol):
#
# Plot the magnitude (in dB) of the transfer function
# to the node numbered '4' in the above netlist
#
  semilogplot(  20*log10(abs(eval(g1, s=I*omega))),
               omega=1..1e09,
               labels=["rad/sec", "Magnitude [dB]"],
               labelfont=[times, bold, 14],
               axes=boxed,
               gridlines,
               labeldirections=[horizontal, vertical]
             );
#
# Plot the phase (in degrees) of the transfer function
# to the node numbered '4' in the above netlist
#
  semilogplot( argument(eval(g1, s=I*omega))*360/2/Pi,
               omega=1..1e09,
               labels=["rad/sec", "Phase [deg]"],
               labelfont=[times, bold, 14],
               axes=boxed,
               gridlines,
               labeldirections=[horizontal, vertical]
             );

#
# Or do the plots the "easy" way by using the BodePlot
# option from the DynamicSystems package - which
# (asmittedly) does a pretty good job of setting more
# or less all the options one might want "correctly"
#
  with(DynamicSystems):
  BodePlot( TransferFunction( g1 ), range=1..1e09);

which produces the following plots

@Fateh 

Your original problem contains

1. the independent variable 'x',
2. the dependent variable 'p(x)' and
3. the parameters 'M', 'k', 'N', and 'lambda'

Your most recent post shows a graph parameterised by a variable 'zeta'.

1. Where did the variable with name 'zeta' come from and how is it related to any of the quantities specified in your original system -ie  'M', 'k', 'N', and 'lambda'?
2. When producing this graph what were the other values of the other parameters in your system?

Maybe you want me to guess what 'zeta' might be because you think I have nothing better to do??

The code below (and attached) shows my best guesses so far, but I won't be making any further contributions based on guesses until you can accurately describe your problem.

The code below (and in the attached odeProb2.mw )

  restart;
  h:= k - (k - 1)*x;
  DE:= diff(p(x), x)
       = -2*M^2
             *
          ( -exp(-M*h/sqrt(N))*exp(M*h/sqrt(N))*sqrt(N)*cosh(M*h/sqrt(N))
            + M*sinh(M*h/sqrt(N))*lambda + exp(-M*h/sqrt(N))*exp(M*h/sqrt(N))*sqrt(N)
          )
          /
          (   exp(-M*h/sqrt(N))
            *
              (   M*exp(2*M*h/sqrt(N))*h 
                - 2*exp(M*h/sqrt(N))*sqrt(N)*cosh(M*h/sqrt(N))
                - M*h + 4*exp(M*h/sqrt(N))*sqrt(N)
                - sqrt(N)*exp(2*M*h/sqrt(N))
                - sqrt(N)
              )
          );
  BC := p(0) = 0, p(1) = 0;
  par:=[0.25, 0.3, 0.4, 0.5, 0.7];
  for j from 1 by 1 to numelems(par) do
      params := [M = 0.2, N = par[j], k =2];
      sol[j] := dsolve(eval({BC, DE}, params), numeric, abserr=1e-03);
      plt[j]:= plots:-odeplot(sol[j], [x, p(x)], x = 0 .. 1);
  od:
  plots:-display( seq(plt[j], j=1..numelems(par)));
  printf( "%12s  %12.8s\n", "N", "lambda");
  seq( printf( "%12.8f  %12.8f\n", par[j], rhs(sol[j](0.5)[3])), j=1..numelems(par));
  

produces the plot (which looks vaguely like your requirement apart from the vertical scale)

and the printout

           N              lambda
  0.25000000    0.65094661
  0.30000000    0.65350011
  0.40000000    0.65672680
  0.50000000    0.65868181
  0.70000000    0.66093406

Note that I have no reason at all to associate your parameter 'zeta' with the parameter 'N'  in your original problem, and also no reason at all to pick the values M=0.2, k=2. I have to make these guesses becuase you are completely incapable of accurately specifying your problem

@Sradharam 

In the last worksheet you posted here  help.mw , you stated

In other words theta(infinify)=1

You now state

"according to our equation(ODE) when xi goes to infinity then the value of theta will be zero"

So which of your two contradictory statements is correct?

This problem is already difficult to solve - when you cannot precisely state the problem, it becomes impossible

why your original code doesn't work, but the attached contains a "workaround" which does work

the code

  restart;
  with(plots):
  interface(version); 

  x1,y1,x2,y2,x3,y3:=0,-3,3,1,5,-2:   
  A := [x1, y1]:
  B := [x2, y2]:
  C := [x3, y3]:

  Barycentre := proc (A, B, t)
                     description "Barycentre de 2 points A(1) et B(t) dans le rapport t";
                     return [(1-t)*A[1]+t*B[1], (1-t)*A[2]+t*B[2]]
                end proc:

  ellip := proc (r1, r2)
                 local a, b, c, d, e, f, D, E, F, eq1, eq2, eq3,
                       eq4, eq5, eq6, x0, y0, EE, r3, sol, Ff, Tg;
                 global A, B, C;
                 r3 := -1/(r2*r1);
                 D := Barycentre(C, B, 1/(1-r1));
                 E := Barycentre(A, C, 1/(1-r2));
                 F := Barycentre(B, A, 1/(1-r3));
                 Ff := proc (x, y)
                             options operator, arrow;
                             a*x^2+2*b*x*y+c*y^2+2*d*x+2*e*y+f
                       end proc;
                 Tg := proc (x0, y0, x, y)
                             options operator, arrow;
                             a*x*x0+b*(x*y0+y*x0)+c*y*y0+d*(x+x0)+e*(y+y0)+f
                       end proc;
                 eq1 := Ff(D[1], D[2]);
                 eq2 := Ff(E[1], E[2]);
                 eq3 := Ff(F[1], F[2]);
                 eq4 := Tg(F[1], F[2], x1, y1);
                 eq5 := Tg(D[1], D[2], x2, y2);
                 eq6 := Tg(E[1], E[2], x3, y3);
                 sol := op(solve([eq1, eq2, eq3, eq4, eq5, eq6], [a, b, c, d, e]));
                 assign(sol);
                 EE := subs(f = 1, Ff(x, y) = 0)
           end proc:

  ellip(-1, -7):
  tri := plot([A, B, C, A], color = blue):
 
  po := plot([A, B, C],
              style = point,
              symbolsize = 15,
              symbol = solidcircle,
              color = red
            ):
  tp := textplot( [ [A[], "A"],
                    [B[], "B"],
                    [C[], "C"]],
                    'align' = {'above', 'left'}
                ):
  ELL := seq
         ( implicitplot
           ( ellip
             ( -7/11-(1/11)*j,
               -1/17-3*j*(1/17)
             ),
             x = 0 .. 5,
             y = -3 .. 1,
             color = ColorTools:-Color
                                 ( [ rand()/10^12,
                                     rand()/10^12,
                                     rand()/10^12
                                   ]
                                 )
          ), j = 1 .. 17):
  display( [tri, ELL, po, tp],
           view = [-.5 .. 5.5, -4 .. 1.5],
           axes = none,
           scaling = constrained,
           size = [500, 500]
         );
  z:=ellip(r1,r2):
  Explore
  ( implicitplot
    ( z,
      x = 0 .. 5,
      y = -3 .. 1
    ),
    parameters = [ r1 = -2.18 .. -.7,
                   r2 = -3 .. -.23]
  );

produces

as well as the "Explore" graph with sliders, which (regrettably) cannot be converted to a graphic form which I can upload here - so you will have o take my work for it, or open the worksheet ellip.mw

I also "tidied up" a few things in this worksheet

@Sradharam 

is the value of the parameter Pr?

In the attached, just as an experiment I produced solutions for Pr=0.1 and "infinity"=10 using Maple's built-in BVP solvers.

Playing around in this worksheet with a few alternatives - higher values for "infinity" and different values for parameter 'Pr' convinced me of two things

1. This new system probably has multiple solutions (just like your earlier one), so the process of setting up the shooting method is probably going to be equally difficult. Perhaps even more so, because the shooting method has to convert two boundary conditions at "infinity" to two initial cnditions D(f)(0) and D(theta)(0)
2. Depending (mainly?) on the value of the parameter 'Pr' then this new system can become numerically unstable (even reported singularies) very quickly!! Apart from anything else, this makes setting up a "shooting" method very difficult, becauseduring the optimization process, the optimiser may try values on some iterations which triggers the numerical instability

Anyhow for what it is worth (and it isn't much!!) you can play around with the attached

odeProb3.mw

@Sradharam 

the DirectSearch packag from here

https://www.maplesoft.com/applications/view.aspx?SID=101333

If I remember correctly, installation is obvious and pretty bulletproof

I probably should have mentioned that the worksheet odeProb2.mw uses a "shooting technique" with an "RKF45 method"

@Preben Alsholm 

and a big upVote.

I spent a lot of time "playing" with your original worksheet and came up with several observations/workarounds

1. I used a value of infinity=10 in my original response - this is waaaay too small
2. In order to achieve a unique value of D(f)(0) which corresponded to f(infinity)=0 (ie the shooting ocndition), for all of the OP's  required parameter combinations (ie  [ alpha=1.5, k=[0, 3, 10, 50] ] and [ alpha=[1, 2, 3, 4], k=1.5] ) then the minimum value for "infinity" was around 200. Some parameter combinations, one can get away with much less, around 30 or so, but to guarantee all of them provide a unique solution, then infinity=~200 seems to be necessary
3. Having established that one needed infinity=~200 to work for all possible combinations, I tried solving it as a BVP with the condition f(infinity=200)=0 for all cases - but generally received the dreaded "Initial Newton iteration is not converging" error message. I tried the usual workarounds for this error message - mainly coming up with possibilities for the 'approxsoln' option, but failed miserably. Not saying it can't be done - just that I couldn't do it. So I decided that this was a limitation of Maple's BVP solvers and looked for an alternative
4. Next possibility was to seriously apply the shooting method - as in convert the BVP to an IVP. The difficulties were
   1. To find the value of D(f)(0) which corresponded to f(infinity=200)=0
   2. To verify that this value was unique, since Preben has already shown that multiple values are possible
5. It is almost impossible to determine whether a unique solution exists using the simple fsolve() command. For anything other that polynomials, fsolve will return one solution, and one has no idea whether other soltuons might exist.  Assuming that one has the DirectSearch package loaded, then a better option is to use its SolveEquations() command, with the AllSolutions=true option. This can in fact be used to generate all three possible solutions for the specific (alpha, k)-combination  which Preben used in the worksheet above. (Obviouly if you don't have the DirectSearch package installed - you are screwed!!)
6. Setting "infinity=200" and having to compute the results of an IVP at f(infinity) results in a maximum-function-evaluation problem, so one has to set maxfun=0 in the dsolve() commands to disable the maxfun limit
7. Applying all of the above I come up with the attached worksheet - even now for certain parameter (ie alpha, k values) I am unconvinced by some of the solutions . The combinations which "bother" me are  [alpha=1.5, k=10], [alpha=2, k=1.5] and [alpha=4, k=1.5] - for reasons which should be "obvious" - the plots of g(xi) just look wrong

Feel free to read/examine the following, but don't even think about re-executing unless

1. You have the DirectSearch package installed
2. You have time to burn - this worksheet takes aroubnd 90minutes to execute on my (admittedly aging) amchine

odeProb2.mw

Forgot the attachment
odeProb19.mw

@Sradharam 

(with the addition of version confirmation) running in Maple2019

@vv 

but if one supplies no labels when creating a dataframe, then default numerical labels will be applied - although these are just row and column numbers. So numerical labels are "allowed" - they have to be. The inconsistency is in their interpretation

@brian bovril 

would be not to use names such as `2` in the data list in the first place - which is (one reason) why Carl's suggestion of structuring the data as a table of records is a good idea.

You can avoid the use of apostrophes in the call to the getter() function in my original code by changing it to

  restart;
#
# data
#
  S := [`206` = Record(mu = 508.001018040,  sigma = 125.002863204708),
          `4` = Record(mu = 1008.001018040, sigma = 167.707232430134),
          `2` = Record(mu = 1208.001018040, sigma = 141.512246146314),
          `5` = Record(mu = 808.001018040,  sigma = 117.156800098735)
       ]:
#
# "getter" function
#
  get:= (data, index, field)-> field = rhs( data[ ListTools:-Search( convert(index, name), lhs~(data)) ] )[field]:
#
# Examples
#
  get(S, 2, mu);
  get(S, 5, sigma);

It looks very much as if your data would be better organised as a table of records, rather than a list.

However ( just to demonstrate that there is always a way! ), you could try the following

  restart;
#
# data
#
  S := [`206` = Record(mu = 508.001018040,  sigma = 125.002863204708),
          `4` = Record(mu = 1008.001018040, sigma = 167.707232430134),
          `2` = Record(mu = 1208.001018040, sigma = 141.512246146314),
          `5` = Record(mu = 808.001018040,  sigma = 117.156800098735)
       ]:
#
# "getter" function
#
  get:= (data, index, field)-> field = rhs( data[ ListTools:-Search( index, lhs~(data)) ] )[field]:
#
# Examples
#
  get(S, `2`, mu);
  get(S, `5`, sigma);

This will return

                      mu = 1208.001018040

                    sigma = 117.156800098735

@vs140580 

click on the plot and use the drag handles??

Unfortunately the 'size' option is not available for 3D plots:-(

Alternatively if you want to plot over a wider range of x, y, z, then you can use the view option, as shown below

with(plots):
pts:= Matrix( [ [ 19.8, 12,    1.62 ],
                [ 24.7, 14,    1.2  ],
                [ 26,   16,    1.25 ],
                [ 35.1, 20,    0.84 ],
                [ 37.8, 20.83, 0.74 ],
                [ 10.8, 9,     3.98 ],
                [ 16.1, 12,    2.46 ],
                [ 38.3, 18,    0.52 ],
                [ 24.5, 16,    1.42 ],
                [ 15.8, 11,    2.35 ],
                [ 28.2, 17.8,  1.19 ],
                [ 56.8, 27,    0.37 ],
                [ 25.4, 17,    1.34 ],
                [ 29.3, 18.1,  1.12 ]
             ]
           ):
  display( [ pointplot3d
             ( pts,
               symbol=solidsphere,
               symbolsize=20,
               color=red
             ),
             pointplot3d
             ( pts,
               style=line,
               color=blue,
               thickness=5
             )
           ],
           labels=["X", "Y", "Z"],
           labelfont=[times, bold, 16],
           view=[0..70, 0..30, 0..5]
         );

@mapleatha 

in Maple18 animMaple18.mw

in Maple 2015 animMaple2015.mw

in Maple 2016 animMaple2016.mw

in Maple 2017 animMaple2017.mw

in Maple 2018 animMaple2018.mw

in Maple 2019 animMaple2019.mw

in Maple 2020 animMaple2020.mw

in MAple 2021 animMaple2021.mw

I'm afraid that I have no conveneinent way of checking the function of this code  in Mpale versions older than Maple 18 (released in 2014). I only have the last eight versions of Maple loaded on this machine. Any version more than eight years old, I can't really help with (other than advise you to update!)