is to use the plottols:-annulus() command, as in

plots:-display( plottools:-annulus( [0,0], 4..5, color=red));

which will produce

precisely what you want or why, but you might want to look at the code below

plots:-display( [ plot( sin(t), t=-Pi..Pi, color=red),
                  plot( 2*sin(t), t=-2*Pi..2*Pi, color=blue),
                  plot( 4*sin(t), t=-4*Pi..4*Pi, color=green)
                ]
              );
plots:-display( [ plot3d( x*sin(t), x=0..1, t=-Pi..Pi,     color=red, style=surface),
                  plot3d( x*sin(t), x=0..2, t=-2*Pi..2*Pi, color=blue, transparency =0.7, style=surface),
                  plot3d( x*sin(t), x=0..8, t=-4*Pi..4*Pi, color=green, transparency =0.7, style=surface)
                ]
              )

which produces the plots (both of which look better in a Maple worksheet than they do on this site)

as far as I can tell the main problem is that you have two nested 'for' loops whichboth  use the same index variable 'i'

for i from 1 to npaths do
      S[0]:=S0 : #` initialise`
      lnS:=ln(S[0]):
      Stot := 0.0:
      for i from 1 to nsteps do
      ........
      od
.......
od

Now every time the 'inner' for loop executes, on exit, the loop index 'i ' will be assigned to nsteps+1 - which in your code will be 2. Unfortunately this will also reset the value of the 'outer' loop index 'i' to 2 and hence the 'outer' loop will never terminate.

The obvious solution is to rename the inner loop index ( 'j' for example!)

I would have attempted this but as far as I can tell nothing in the 'outer' loop depends on the loop index 'i' - whihc presumably means that something in the 'inner' loop depends on the 'inner' loop index and sometihng in the 'loop' depends on the outer loop index. No way I can tell which is which.

Just change the above code to

for i from 1 to npaths do
      S[0]:=S0 : #` initialise`
      lnS:=ln(S[0]):
      Stot := 0.0:
      for j from 1 to nsteps do
      ........
      od
.......
od

and rewrite the dependencies in the 'inner' loop appropriately - ie what depends on the 'outer' loop variable 'i' and what depends on the 'inner' loop variable 'j'



 

since all you are doing is setting up two names for the same numerical value - ie in your example

A[1] will return 14.9057064333276, and

Y[1,1] will return 14.9057064333276

so why bother it appear to be a complete waste of computation time and memory!!

Can it be done? Of course. Two ways are shown below (and then are probably others

  restart;
  A:=Vector[row](19, [ 14.9057064333276, 14.4384716751962, 14.0155569170648, 13.6381346589334,
                       13.3075724008020, 13.0254476426706, 12.7935628845392, 12.6139606264079,
                       12.4889383682765, 12.4210636101451, 12.4131888520137, 12.4684665938823,
                       12.5903643357509, 12.7826795776196, 13.0495548194882, 13.3954925613568,
                       13.8253703032254, 14.3444555450940, 14.9584207869626]
                );
#
# Example
#
  A[12];
#
# One way, with example
#
  seq( assign( Y[1,j], A[j]), j=1..numelems(A)):
  Y[1,12];
#
# Another way, with example
#
  Z:=convert(A, Matrix):
  Z[1,12];

whihc will produce

                        12.4684665938823

                        12.4684665938823

                        12.4684665938823

You should try to understand the difference between "local" and "global" variables, as well as "scoping" rules (ie when a variable is considered "local" to a function and when it is considered "global").

If you cannot grasp this concept, you may find it helpful to rewrite code to avoid "apparent" name clashes between local and global variables. For example the following

  restart;
  f := aLocalVariableName1 -> piecewise(0 <= t, 1, t < 0, t - 1);
  F := aLocalVariableName2 -> int(f(anyDummyVariableName), anyDummyVariableName = 1 .. aLocalVariableName2);
  F(x);

will produce

piecewise( t < 0, t*(x - 1) + 1 - x,
                  0 <= t, x - 1)

whihc is (probably?) what you expect. But you might have to think about the difference between this code and your original

"numerical" approaches to this problem are still being proposed.

In an earlier answer I pointed out that the problem can (I admit somewhat surprisingly)  be solved symbolically. This is so much more efficient than the numerical approaches.

For the numerical approach proposed by CarlLove, (s)he reports the cost for the various stages as

memory used=8.47GiB, alloc change=0 bytes,
cpu time=84.97s, real time=78.92s, gc time=14.45s

memory used=4.69GiB, alloc change=48.00MiB, 
cpu time=50.72s, real time=47.10s, gc time=8.64s

For the numerical approach  by Acer, (s)he  reports the cost as

memory used=8.44GiB, alloc change=36.00MiB, cpu time=64.89s,
real time=64.91s, gc time=5.06s

For the symbolic approach which I advocate the cost is

memory used=161.75MiB, alloc change=136.00MiB, cpu time=1.92s, real time=1.87s, gc time=140.40ms
memory used=0.77MiB, alloc change=0 bytes, cpu time=140.00ms, real time=139.00ms, gc time=0ns
memory used=42.19MiB, alloc change=0 bytes, cpu time=515.00ms, real time=489.00ms, gc time=93.60ms

So the symbolic appraoch is ~200X more efficient on memory usage and ~30X more efficient on cpu time.
Would someone please explain very clearly why the advocate a "numerical" approach when the "symbolic approach is so much more efficient - because I'm not getting it!

The code shown below, which produces the timings/memoryUsage shown above

  restart:
  with(plots):
#
# Define system
#
  local delta:
  assume(delta,real):
  params:=[phi=0, lambda=0.1, N=5, M=sqrt(N*(N+1))*exp(I*phi), omegap=10]:
  dsys:={ diff(n(t),t)=-2*(n(t)-N)+(u(t)-M)*exp(-2*I*omegap*t/lambda)+((z(t)-conjugate(M))*exp(2*I*omegap*t/lambda)),
          diff(u(t),t)=-2*(1-I*delta)*u(t)+2*(n(t)-N)*exp(2*I*omegap*t/lambda)+2*M,
          diff((z(t),t))=-2*(1+I*delta)*z(t)+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)+2*conjugate(M)}:
#
# Solve the system symbolically
#
  ans2:=CodeTools:-Usage(dsolve(eval[recurse]( dsys, params) union {n(0)=0,u(0)=0,z(0)=0})):
#
# Produce a 3D-plot as a function of t and delta from
# the analytic solution
#
  CodeTools:-Usage
             ( plot3d
               ( eval(n(t), ans2),
                 t=-1..1,
                 delta=-10..10,
                 style=surface
               )
             );
  CodeTools:-Usage
             ( plot
               ( eval( eval(n(t), ans2),t=0.5),
                 delta=-10..10,
                 style=surface
               )
             );

will produce the graphs below (shown by previous posters) - just a lot faster and more efficiently

Worksheet for the above is below

odeAnal2.mw

cosider the code

 f__1 := sqrt(4*a^2 + lambda__g^2)*c/(2*lambda__g*a);
  f__2 := c*sqrt(1/lambda__g^2 + 1/(4*a^2));
  is(f__1=f__2);
  is(f__1=f__2) assuming a>0, lambda__g>0;

where the two 'is()' commands, return false and true respectively

(at least to me!), this system can be solved analytically

The code below

1. Sets delta=1, solves the system analytically and plots n(t)
2. Solves the system with delta *unknown*, then plots the solution for delta=1 (mainly for comparison with the above) and also produces a 3d plot of n(t) versus 't' and 'delta'

A couple of points

1. The analytic solution is pretty ugly - but it exists
2. The solution for n(t) appears to be always real, but that for z(t) is complex. I haven't tried to produce any plots for z(t)

  restart:
  with(plots):
#
# Define system
#
  local delta:
  assume(delta,real):
  params:=[phi=0, lambda=0.1, N=5, M=sqrt(N*(N+1))*exp(I*phi), omegap=10]:
  dsys:={ diff(n(t),t)=-2*(n(t)-N)+(u(t)-M)*exp(-2*I*omegap*t/lambda)+((z(t)-conjugate(M))*exp(2*I*omegap*t/lambda)),
          diff(u(t),t)=-2*(1-I*delta)*u(t)+2*(n(t)-N)*exp(2*I*omegap*t/lambda)+2*M,
          diff((z(t),t))=-2*(1+I*delta)*z(t)+2*(n(t)-N)*exp(-2*I*omegap*t/lambda)+2*conjugate(M)}:
 

#
# Solve the system analytically for the case delta=1 (just as
# an example)
#
  ans1:=dsolve(eval[recurse](dsys, [ params[], delta=1] ) union {n(0)=0,u(0)=0,z(0)=0}):
  p1:= plot( eval(n(t),ans1),
             t=0..10,
             color=red
           );
#
# Solve the system analytically for delta "unknown"
#
  ans2:=dsolve(eval[recurse]( dsys, params) union {n(0)=0,u(0)=0,z(0)=0}):
#
# Evaluate this solution at delta=1 and plot - just for
# comparison with the above figure
#
  plot(  eval(eval(n(t), ans2), delta=1),
         t=0..10,
         color=red
      );
#
# Produce a 3D-plot as a function of t and delta from
# the analytic solution
#
   plot3d( eval(n(t), ans2),
           t=0..10,
           delta=0..10,
           style=surface
         );

produces the plots

More details in

odeAnal.mw

For the data you supply, the code below shows that the fitting parameters are

[X__0 = 127.116936525446, Ym = 0.668691089920431, c = 121.021223406666, d = -0.134524707088629]

and the condition Y=Ym/2 is met when X is one of [36.39639448, 679.0013662]

Basic code below, executable worksheet at the end

  restart;
  with(Statistics):
  with(RootFinding):
#
# Define the function and the data to fit
#
  func:=Y = Ym*(c/((d + 1)*(X - X__0) + c))^(d + 1)*exp((X - X__0)*(d + 1)^2/((d + 1)*(X - X__0) + c)):
  data:= < <20   | 0.188081>,
           <30   | 0.193294>,
           <60   | 0.561155>,
           <110  | 0.688284>,
           <160  | 0.607743>,
           <410  | 0.478664>,
           <710  | 0.35032>,
           <1010 | 0.213999>,
           <1410 | 0.204433>
         >:

#
# Fit the data to the function. Note that for functions this
# complicated it is a very, VERY good idea to supply some
# sort of initial "guess
#
  myFit:= Fit( rhs(func),
               data,
               X,
               initialvalues=[Ym = 1/2, X__0=100, c = 100],
               output=parametervalues
             );
#
# Plot the input data and the fitted function
#
  plots:-display( [ plot( eval(rhs(func), myFit), X=0..1500),
                    plot( data[..,1], data[..,2], style=point)
                  ]
                 );
#
# Evaluate the original function with the fitted parameters
# just as some kind of "sanity check"
#
  eval[recurse](func, [ myFit[], Y=Ym/2 ]);

#
# Convert the fitted function to a function with "roots" by
# constructing rhs(fit)-lhs(fit) - now we are interested in
# zero-crossings
# 
  f:=unapply(rhs(%)-lhs(%), X):
#
# This plot shows that there appear to be two "roots
#
  plot(f(X), X=0..1500);
#
# Find both of these zero-crossings
#
  rt1:=0:
  rts:=Array(1..2):
  for i from 1 by 1 to 2 do
      rts[i]:=NextZero(f, rt1, maxdistance=1000);
      rt1:=rts[i];
  od:
#
# Display the required roots, ie the values of 'X' at
# which Y=Ym/2 in the original function
#
  rts;
#
# A check! Display the value of the parameter Ym
#
  eval(Ym, myFit);
#
# Evaluate the original function with the fitted
# parameters at the values of the variable 'X'
# obtained as roots. These should be 1/2 of the
# value of Ym obtained above
#
  eval( func, [myFit[], X=rts[1]]);
  eval( func, [myFit[], X=rts[2]]);

rtFind.mw

Running a couple of simple tests with different values for the "parameters" I can come up with cases where there are two solutions and cases where there are three solutions (and I'm not claiming that these are the only cases!).

The code below

  restart;
  with(RealDomain):
#
# Equation of interest
#
  eqn1:=Y = Ym*(c/((d + 1)*(X - X__0) + c))^(d + 1)*exp((X - X__0)*(d + 1)^2/((d + 1)*(X - X__0) + c));
#
# Evaluate equation for Y=Ym/2;
#
  eqn2:= eval(eqn1, Y=Ym/2):
#
# Try an "arbitrary" set of values for c, d, X__0. ( Exact
#  knowledge of these parameter values would be *very* useful!)
#
  params:=[c=1, d=1, X__0=1]:
#
# Hmm 'solve()' produces three(!) solutions - maybe a negative
# value for 'X' is not meaningful?
#
  sol:=evalf~
       ( [ solve
           ( eval
             ( eqn2,
               params
             ),
             X
           )
         ]
       );
#
# Check the value of 'Y' for each of these three solutions
#
  seq( eval(eqn1, [params[], X=sol[j]]), j=1..3);
#
# A different set of "arbitrary" values for these parameters
# produces a solution dependent on a binary-valued parameter
#
  params:=[c=1, d=2, X__0=3]:
  sol:= evalf~
        ( [ solve
            ( eval
              ( eqn2,
                params
              ),
              X
            )
          ]
        );
#
# So evaluate at both binary values (ie 0,1)
#
  vals:=[seq( eval(sol, indets(sol, name)[]=j)[], j in [0,1])];
#
# Check the value of 'Y' for each of these two solutions
#
  seq( eval(eqn1, [params[], X=vals[j]]), j=1..2);

will generate and verify the soutions for two different lists of parameter values.

1. The list [c=1, d=1, X__0=1], will produce the solutions [-1.872951360, 1.813645761, 0.7406205565], for which the verification step produces Y = 0.4999999995*Ym, Y = 0.4999999998*Ym, Y = 0.5000000001*Ym
2. The list [c=1, d=2, X__0=3] will produce the solutions [3.383956320, 2.847650669], for which the verification step produces  Y = 0.4999999991*Ym, Y = 0.4999999940*Ym

 so all appear to be "correct" solutions for Y=Ym/2. See also the attached

getSol.mw

The code below

  restart;
  with(ExcelTools):
  fname:="C:/Users/TomLeslie/Desktop/":
  L:=[[[1,2],3],[[1,3],4],[[1,4],4],[[1,7],7]];
  M:=Matrix( [ ["A", "B"],
               seq( [convert(L[j][1],string), L[j][2]], j=1..numelems(L))
             ]
           );
  Export( M,
          cat(fname, "dat.xlsx"),
           "Sheet1",
          "B2"
        );

will produce the Excel file attached

dat.xlsx

You will have to change the string "fname" to something appropriate for your machine

the code below?

  restart;
#
# Define function
#
  Crunch:= (x::numeric)-> `if`(x>=100, x/100, x+Crunch(10*x));
#
# Apply function to (four) list entries
#
  Crunch~([ 2, 202, 200, 2000]);

whihc willl produce

[24, 101/50, 2, 20]

which isn't particularly elegant!

  restart;
#
# Generate a random array
#
  r:=rand(1..10):
  A:=Array( [seq(r(), j=1..20)]);
#
# indices to remove
#
  inds:={3, 7, 8, 12, 15};
#
# Utility function
#
  f:=( x::Array, y::set ) -> local i,j:
                             [ seq( x[j],
                                    j in {$i=1..upperbound(x)} minus y
                                  )
                             ];
#
# New Array
#
  B:= Array( f(A, inds) );

rmArr.mw

Increasing Digits to 15 makes the problem "go away"

The code

restart;
Digits:=15:
with(plots):
with(Statistics):
X1 := RandomVariable(BetaDistribution(4, 4));
X2 := RandomVariable(BetaDistribution(4, 4));
plot([PDF(X1, t), PDF(0.5*X1 + 0.5*X2, t)],
      t = 0 .. 1,
      color = [red, blue],
      legend = ["PDF X1", "PDF average"]
    );

produces

as in

esm := ExponentialSmoothingModel(ts);
pList := GetParameters(esm);
eval(alpha, pList);
eval(beta, pList);
GetParameter(esm, [alpha, beta]);

See the worksheet here expSmooth.mw for more detail