You mentioned it, but this is not how your code was written. You started j from 2 and the code uses values at time j to compute values at time j+1. As my code showed, you need to start with j=1. If you insist on starting with j=2, then you need to decrease all second arguments by 1 (use data at time j-1 to compute values at time j).

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

Why do you say this "does not work at all"?

Is it because these loops produce no output? This is caused by two facts. First, the colon in the second end do suppresses all output. But, even changing this to a semi-colon does not produce any output. This is because Maple does not, by default, display output produced at too deep of a level. To change this, execute the command:

printlevel := 10;

See also the online help, where you will see how the "level" is determined and that using printlevel=1000 is not unreasonable for many development situations.

You say that f[n,m] has been setup. All you really need is the information at the initial timestep. You start j at 2, and compute the solution at time j+1. I'm guessing you really want to start j at 1 (or 0).

Here's what I did to see that your code is not unreasonable:

f[1,1],f[2,1],f[3,1],f[4,1]:=0,1,2,3; # fictitious data
m,n,d := 3,4,0.1;                     # just for testing
for j from 1 to m-1 do
       for i from 1 to n do
 if i=1 then
       f[i,j+1]:=f[i,j]+d*(f[i+1,j]-2*f[i,j]+f[n,j])
 elif i=n then
       f[i,j+1]:=f[i,j]+d*(f[1,j]-2*f[i,j]+f[i-1,j])
 else
       f[i,j+1]:=f[i,j]+d*(f[i+1,j]-2*f[i,j]+f[i-1,j]);
 end if;
 end do; end do;

I would probably rewrite the inner loop without the conditional, as follows

for j from 1 to m-1 do
  f[1,j+1]:=f[1,j]+d*(f[2,j]-2*f[1,j]+f[n,j]);
  for i from 2 to n-1 do
    f[i,j+1]:=f[i,j]+d*(f[i+1,j]-2*f[i,j]+f[i-1,j]);
  end do;
  f[n,j+1]:=f[n,j]+d*(f[1,j]-2*f[n,j]+f[n-1,j])
end do;

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

You have to read the help carefully, paying close attention to the syntax and assumptions that are made.

In this case, we are talking about the Window "element". Every window has to have a layout, you can do this by saying 'layout'=, but if you omit this label, then any list is automatically assumed to be the layout. So, when you entered (note that there was an extra ' in the Button element)

with( Maplets[Elements] ):
DoubleX:=Maplet(Window('title'="DoubleX",
  ["Enter X",TextField['TF']()],
  TextBox['XDoubled']('editable'='true'),
  Button("XDoubled",Evaluate('XDoubled'=2*'TF'))
)):

the layout was taken to be ["Enter X",TextField['TF']()], which creates an (unnamed) label and the TF textbox in a column. The XDoubled textbox and the button are not recognized as part of the layout, but do not generate any error message.

Now, when you entered (your code had a missing comma in addition to the extra ')

DoubleX:=Maplet(Window('title'="DoubleX",[
  ["Enter X",TextField['TF']()],
  TextBox['XDoubled']('editable'='true'),
  Button("XDoubled",Evaluate('XDoubled'=2*'TF'))
] )):

now the window consists of two items, the title and the layout (everything inside the outermost square brackets). That is the difference.

To be safe, at least to start, you should write out every option to each element. (I wish the examples in the help system would do a better job of this). The help is complete and carefully explains how each implicit option is handled, but this requires the careful reading that I mentioned at the start.

Here is an excerpt from the Maplets,Elements,Window help page with the relevant information for this example:

The Window element features can be modified by using options. To simplify specifying options in the Maplets package, certain options and contents can be set without using an equation. The following table lists elements, symbols, and types (in the left column) and the corresponding option or content (in the right column) to which inputs of this type are, by default, assigned.

Elements, Symbols, or Types                       Assumed Option or Content
                                                                           
list, BoxLayout element, or GridLayout element    layout option            
MenuBar element                                                        menubar option           
refID                                                                               reference option         
string or symbol                                                           title option             
ToolBar element                                                          toolbar option

Note, in particular, that you could omit the 'title'= since a string is automatically assumed to be the title.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

You wrote:

Many thanks for this Doug. Your post has helped and I have read the help but I'm still missing something. I have a list of the following points:

[[3, 3], [5, 20], [7, 63], [9, 144], [11, 275]]

When I use interactive as you describe, I get the options as you mention. Using these I can get a simple curve through all these points which  I am happy with. However, when I then try to plot this, the resulting curve looks completely different. Trying something similar but with Spline instead of interactive just resulted in a straight line plot (y=x)

What do you mean "when I then try to plot this"? What, exactly, did you try to plot? How did it look different?

What did you do when you tried "something similar but with Spline"? What commands did you use? Did the curve look linear or did you get an equation? Was it a "straight line" or was it exactly "y=x"?

The answers to these questions will help us to know exactly where to direct our efforts. It would probably be easiest if you uploaded your worksheet.

Here is something I did with your data:

with( CurveFitting ):
data := [[3, 3], [5, 20], [7, 63], [9, 144], [11, 275]];
               [[3, 3], [5, 20], [7, 63], [9, 144], [11, 275]]
f := Interactive( data ); # choose the fit you want
                              1  3   1  2   1  
                              - x  - - x  + - x
                              4      2      4  
g := Interactive( data ); # choose another fit
                               9664     34805  2
                             - ---- x + ----- x 
                               685      10001   
plot( [data,f,g], x=-2..15,
      style=[point,line,line], color=[red,blue,green],
      legend=["data pts","polynomial interp","least squares parabola"] );

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

I'm glad this was helpful. I meant to make explicit mention to the fact that I changed the time interval in each DEplot command. If I left the time interval starting at t=0, then Maple would draw the solutions backward in time from t_0 to t=0 (and forward from t_0 to t=10. If this is what you want, it's no problem to make this refinement.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

All you need to do is change the value of the independent variable in your initial conditions. Here I do 4 different sets of initial conditions, then display the results in a single plot.

restart;
with( DEtools ):
de := diff(y(t), t) = y(t)-y(t)^2+.3*sin(2*Pi*t);
inits0:=[[y(0)=0.5],[y(0)=1.0],[y(0)=1.5]]:
p0 := DEplot(de, y(t), t = 0 .. 10, inits0, y(t) = 0 .. 2, stepsize = .1):
inits1:=[[y(1)=0.5],[y(1)=1.0],[y(1)=1.5]]:
p1 := DEplot(de, y(t), t = 1 .. 10, inits1, y(t) = 0 .. 2, stepsize = .1, linecolor=cyan, arrows=none):
inits1:=[[y(2)=0.5],[y(2)=1.0],[y(2)=1.5]]:
p2 := DEplot(de, y(t), t = 2 .. 10, inits1, y(t) = 0 .. 2, stepsize = .1, linecolor=green, arrows=none):
inits1:=[[y(5)=0.5],[y(5)=1.0],[y(5)=1.5]]:
p5 := DEplot(de, y(t), t = 5 .. 10, inits1, y(t) = 0 .. 2, stepsize = .1, linecolor=brown,arrows=none):
plots[display]( p0,p1,p2,p5 );

Is this useful?

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

Take a look at ?Maplets

There are lots of links from this page to additional information. Most are pretty good, with lots of examples. I would not spend any time looking at the Maplet Builder (see previous comments on MaplePrimes).

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

You can use the while option as I suggested. Also, the previous suggestions to use break are appropriate here. Since you have not provided full information, I cannot fully test what I am about to suggest, but these give the general idea. First, with the break command:

for i from 0 to 2 do
       a[i+1]:=subs(z1=x[i+1][1],z2=x[i+1][2],z3=x[i+1][3],G):
       if a[i+1]=0 then break end if
end do:

And now with the while option:

for i from 0 to 2 while a[i]<>0 do
       a[i+1]:=subs(z1=x[i+1][1],z2=x[i+1][2],z3=x[i+1][3],G)
end do:

The latter one might need some tweaking to get it started (make sure a[0]<>0), but these give the general ideas.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

I think you will want to look at using a module to implement your package. Then, the ModuleLoad procedure within the module is executed when the module (package) is loaded. See the online help pages: ?module and ?ModuleLoad

Does this help?

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

It appears you are trying to find the first positive value where the function f is increasing. There are lots of problems with what you have tried. I will attempt to put you back on the right track.

First, your definition of f does not appear to depend on b. You don't tell us the values you are using for z[1], z[2], and z[3], maybe these involve b. To define a function of the variable b, I would suggest using the arrow notation:

f := b -> (b-4)^4+(b-3)^2+4*(b+5)^4;

To construct a loop that terminates the first time it detects an increase in the value of f, you could do the following:

for b from 0 to 0.003 by 0.0001 while f(b)<f(b-0.0001) do
end do;
if b<> 0.003 then
  rhs_interval := b;
end if;
print( b, f(b), f(b-0.0001), f(b)-f(b-0.0001) );

This loop does nothing except increase the value of b and check that the value of f is smaller than at the previous value of b. It there is a tie, or if the value of f is larger than the previous value, then the loop terminates. When the loop terminates, check if it simply ran out of values to check or if it actually stopped due to the stopping condition in the while clause.

This is not the way I would approach this type of problem. You should be able to use calculus to reformulate the problem in a way that you could ask Maple to solve some collection of equations.

It would help to know more about what you are really trying to do, and what level of mathematics you know and can use for to answer this question.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

The order of the results from fsolve can change. You cannot do anything about this. You should be thinking about a way to work with these results to put them in the format you desire.

The following shows how you can create an ordered pair [u2,fout] from the contents of ddmax:

ddmax := {fout = 0.3873538368e-1, u2 = 9.708916648};
ddmax2 := eval( [u2,fout], ddmax );

If you want them in a different order, or any other expressions involving fout and/or u2, simply put this in the first argument of the eval command, e.g.,

eval( [u2/fout, u2^2+fout^2], ddmax );

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

You don't give us very much definite information about your problem, but let's see if some of the ideas I'll describe are helpful for you.

If you have a collection of data points (in 2D), here's one way to visualize the data and to see some of the different curves based on this data (polynomial interpolant, splines, least squares, ...). All of this is done through an interactive assistant:

with(CurveFitting):
pts := [[1.3, .56], [3.5, 1.1], [6.7, 2.97], [7.9, 4.3], [9.3, 6.1]];
Interactive( pts );

If you don't have the data points, this will call up an interface for entering the data before giving you the different options for fitting the data.

Interactive();

A worksheet with explicit examples and instructions for using the CurveFitting package can be accessed by entering:

?CurveFitting,examples

If none of this leads you to the results you are seeking, please let us know what you are trying to do and give us an idea of what you have tried and how it has not been what you need.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

While your definition of ode might look correct, the RHS is not the function you expect it to be. Maple understands square brackets as delimiters for a list, not an interval. So, when you wrote t in [0,Pi], this will be true only when t=0 or t=Pi. The simplest fix is to replace t in [0,Pi] with t<Pi

ICS  := y(0)=0, D(y)(0)=0:
ODE2 := diff(y(t),t$2)+y(t)=piecewise( t<Pi, t, t>Pi, Pi );
         / d  / d      \\                                          
         |--- |--- y(t)|| + y(t) = piecewise(t < Pi, t, Pi < t, Pi)
         \ dt \ dt     //                                          
dsolve( {ODE2,ICS}, y(t) );
       y(t) = piecewise(t < Pi, -sin(t) + t, Pi <= t, -2 sin(t) + Pi)

To see the difference in the right-hand sides, consider

plot( [rhs(ODE),rhs(ODE2)], t=-1..5, style=point, symbol=[circle,cross] );

The value of the RHS for t<0 is not important to this problem. But, if you want to insist that the RHS is zero (e.g.) for t<0, then insert  t<0, 0 before the current first argument in the piecewise command.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

The only thing you are doing "wrong" is that you did not include the arguments of exp_x in the plot3d call. This line should be:

plot3d(exp_x(y,z),y=-15..15, z=-15..15);

Now, the plot looks pretty rough because outside the circle with radius 5, exp_x(y,z) is an imaginary number that cannot be plotted (by plot3d). You have a couple of options for improving the plot. First, you could reduce the domain to

plot3d(exp_x(y,z),y=-5..5, z=-5..5);

This plot still is rather jagged at the edges. This is becuase Maple uses a rectangular grid to sample points for the plot. You can improve the resolution by adding the grid= option, as in the following. (I also add the axes=boxed option in antiicipation for the next

plot3d(exp_x(y,z),y=-5..5, z=-5..5, grid=[95,95], axes=boxed );

This looks better, but is still not identifiable as a sphere because of the different scales used on the vertical axis. To force all three axes to be shown with the same scale, add the scalinig=constrained option.

plot3d(exp_x(y,z),y=-5..5, z=-5..5, grid=[95,95], axes=boxed, scaling=constrained );

This should be looking better. Now, let's add the bottom half of the sphere.

plot3d([exp_x(y,z),-exp_x(y,z)], y=-5..5, z=-5..5, grid=[95,95], axes=boxed, scaling=constrained );

The scaling= option is unneeded in this command. (If you stick with the initial intervals of -15..15, then you will need to use grid=[285,285] to get the comparable resolution, and scaling=constrained will be needed.)

I hope this gets you closer to where you need to be.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed

I have now created an initial implementation of routines that plot the little direction arrows in the phase portrait for a system of ODEs. The PParrow1 procedure plots the single arrow in a coordinate dirction and PParrow2 produces the pair of direction arrows at a specified point.

PParrow2 := proc( sys, pt )
  local opt;
  if  _npassed>2 then opt := _passed[3..-1] else opt := NULL end if;
  plots:-display( [ PParrow1(sys[1],pt,<1,0>,opt),
                    PParrow1(sys[2],pt,<0,1>,opt) ] )
end proc:
PParrow1 := proc( eqn, pt, dir )
  local d, p, opt;
  if  _npassed>3 then opt := _passed[4..-1] else opt := NULL end if;
  d := signum( eval(eqn,pt) )*dir;
  p  := map(rhs,pt);
  plots:-arrow(p,d,opt)
end proc: 

To work on this problem, let's enter the complete system of 3 ODEs

SYS := [
  diff(u(t),t) = (B/alpha-B+B*u(t)-phi(t))*u(t),
  diff(phi(t),t) = (1/theta)*((alpha*z(t)-delta+rho)-(z(t)-delta-phi(t)))*phi(t),
  diff(z(t),t) = (alpha-1)*(z(t)-B/alpha)*z(t)
];

and the parameter values provided by the orignal poster:

param := [
  B=1,
  alpha=0.3,
  delta=0.03,
  rho=0.02,
  theta=1,
  z(t)=1/0.3,
  NULL ];

Let's first find the specific system of interest, then find algebraic descriptions for the nullcline curves (implicitly)

SYS2 := eval(SYS[1..2],param);
SYS3 := eval(SYS2,[u(t)=u,phi(t)=phi]);

The nullclines in the u (blue) and phi (green) phase space can be plotted as follows:

nullcline1 := implicitplot( SYS3[1], u=-3..3, phi=-3..6,
                            grid=[65,65], thickness=2, color=blue, legend="u ' = 0" ):
nullcline2 := implicitplot( SYS3[2], u=-3..3, phi=-3..6,
                            thickness=2, color=green, legend="phi ' = 0"  ):
ncP := display( [nullcline1, nullcline2] ):
ncP;

Now, the PParrow2 command will be used to insert the direction arrows at one point in each of the separate regions of phase space created by the nullclines.

display( [ncP,
          PParrow2(map(rhs,SYS3),[u= 1.5,phi= -2],length=0.25),
          PParrow2(map(rhs,SYS3),[u= 1.5,phi=  1],length=0.25),
          PParrow2(map(rhs,SYS3),[u= 1.5,phi=  3],length=0.25),
          PParrow2(map(rhs,SYS3),[u= 1.5,phi=  5],length=0.25),
          PParrow2(map(rhs,SYS3),[u=-1,  phi= -2],length=0.25),
          PParrow2(map(rhs,SYS3),[u=-1,  phi=0.5],length=0.25),
          PParrow2(map(rhs,SYS3),[u=-1,  phi=  2],length=0.25),
          PParrow2(map(rhs,SYS3),[u=-1,  phi=  4],length=0.25)
         ] );

Note that I have assumed z is at its equilibrium solution (z=B/alpha).

You should be able to adapt this to your specific problem, if this is not it.

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.ed