With y1, ..., y4 as you had them:

> series((2*y1-5*y2+4*y3-y4)/(h^2), h);

`@@`(D,2)(y)(0)-11/12*`@@`(D,4)(y)(0)*h^2+O(h^3)

This might work. It is based on the fact that plot(..., discont=true) produces an object of the form CURVES(L1, L2, ..., options) where each Lj is a list of points. The vertical asymptotes should be at x values between those of the last entry in one list and the first in the next. The syntax for vplot is the same as for plot.

> vplot:= proc()
    local P, C, segs, i, p1, p2, x;
    P:= plot(args, discont=true):
    C:= op(1,P):
    if not type(C,specfunc(anything,CURVES)) then 
      return P 
    end if;
    segs:= NULL;
    for i from 1 to nops(C)-1 do
      if type(op(i,C),listlist) and type(op(i+1,C),listlist) 
      then
        p1:= op(i,C)[-1]; 
        p2:= op(i+1,C)[1];
        x:= (p1[1]+p2[1])/2;
        segs:= segs, [[x,p1[2]],[x,p2[2]]];
      end if
    end do;
    plots[display](P, plot([segs],linestyle=3,colour=black));
  end proc;

For example:

> vplot(tan(x),x=0..3*Pi,-10..10);
  vplot(arctan(tan(x)),x=0..3*Pi);

As the error message says, GenerateMatrix expects a list but received a Vector. You could try

> MV:= GenerateMatrix(convert(y,list),convert(x,list));

As you noticed, the result (which I assigned to MV) is a Matrix and a Vector. The Matrix is MV[1].

As the error message says, GenerateMatrix expects a list but received a Vector. You could try

> MV:= GenerateMatrix(convert(y,list),convert(x,list));

As you noticed, the result (which I assigned to MV) is a Matrix and a Vector. The Matrix is MV[1].

For (3): I assume your matrix is hermitian (e.g. real and symmetric), so the eigenvalues should of course be real. However, if it doesn't have shape = hermitian or shape = symmetric, Maple doesn't know this, and the eigenvalues produced numerically by Eigenvalues may sometimes have small nonzero imaginary parts. For example:

> with(LinearAlgebra):
  B:= RandomMatrix(20,20):
  Di := DiagonalMatrix([0$10,1$10]):
  A:= B . Di . Transpose(B):
  convert(Eigenvalues(evalf(A)),list);

[130484.415156932635+0.*I, 122280.763414717862+0.*I, 87936.0963894075394+0.*I, 80208.9073262921884+0.*I, 66275.9901489112817+0.*I, 47562.3730272897446+0.*I, 39306.6050666601441+0.*I, 23571.0830624543960+0.*I, 12502.9521446575400+0.*I, 11631.8142626769386+0.*I, .144642265187316840e-10+.242049577586753941e-12*I, .144642265187316840e-10-.242049577586753941e-12*I, -.609594215585161176e-11+.137938928318574424e-11*I, -.609594215585161176e-11-.137938928318574424e-11*I, .539962298366958798e-11+0.*I, .316339746841558758e-11+.325321480613692524e-11*I, .316339746841558758e-11-.325321480613692524e-11*I, -.962503680057553204e-12+.196442958624270660e-11*I, -.962503680057553204e-12-.196442958624270660e-11*I, -.131265631539004777e-11+0.*I] This problem of imaginary parts can be combatted by something like:

> A1 := Matrix(A, shape=symmetric);
  Eigenvalues(evalf(A1));

[-.300538544341638301e-10, -.159312544468454100e-10, -.116997883938897732e-10, -.377349888401283454e-11, -.974198084677210350e-12, -.592861334731303932e-12, .165311614100773302e-11, .819101484874475244e-11, .114713694168626624e-10, .146458395266595254e-10, 11631.8142626769258, 12502.9521446575128, 23571.0830624544288, 39306.6050666601950, 47562.3730272897010, 66275.9901489112090, 80208.9073262920720, 87936.0963894073648, 122280.763414717993, 130484.415156932562] However, there's still a problem of roundoff error causing the eigenvalues that should be exactly 0 to become nonzero. If (as in this case) the matrix entries are exact rationals, and the matrix is small enough that computing the characteristic polynomial is feasible, you can use a Sturm sequence to avoid trying to compute eigenvalues altogether:

> C:= CharacteristicPolynomial(A, t):
  positives:= sturm(sturmseq(C,t),t,0,infinity);
  zeros:= ldegree(C,t);
  negatives:= degree(C,t)-positives-zeros;
  signature:= positives - negatives;

For (3): I assume your matrix is hermitian (e.g. real and symmetric), so the eigenvalues should of course be real. However, if it doesn't have shape = hermitian or shape = symmetric, Maple doesn't know this, and the eigenvalues produced numerically by Eigenvalues may sometimes have small nonzero imaginary parts. For example:

> with(LinearAlgebra):
  B:= RandomMatrix(20,20):
  Di := DiagonalMatrix([0$10,1$10]):
  A:= B . Di . Transpose(B):
  convert(Eigenvalues(evalf(A)),list);

[130484.415156932635+0.*I, 122280.763414717862+0.*I, 87936.0963894075394+0.*I, 80208.9073262921884+0.*I, 66275.9901489112817+0.*I, 47562.3730272897446+0.*I, 39306.6050666601441+0.*I, 23571.0830624543960+0.*I, 12502.9521446575400+0.*I, 11631.8142626769386+0.*I, .144642265187316840e-10+.242049577586753941e-12*I, .144642265187316840e-10-.242049577586753941e-12*I, -.609594215585161176e-11+.137938928318574424e-11*I, -.609594215585161176e-11-.137938928318574424e-11*I, .539962298366958798e-11+0.*I, .316339746841558758e-11+.325321480613692524e-11*I, .316339746841558758e-11-.325321480613692524e-11*I, -.962503680057553204e-12+.196442958624270660e-11*I, -.962503680057553204e-12-.196442958624270660e-11*I, -.131265631539004777e-11+0.*I] This problem of imaginary parts can be combatted by something like:

> A1 := Matrix(A, shape=symmetric);
  Eigenvalues(evalf(A1));

[-.300538544341638301e-10, -.159312544468454100e-10, -.116997883938897732e-10, -.377349888401283454e-11, -.974198084677210350e-12, -.592861334731303932e-12, .165311614100773302e-11, .819101484874475244e-11, .114713694168626624e-10, .146458395266595254e-10, 11631.8142626769258, 12502.9521446575128, 23571.0830624544288, 39306.6050666601950, 47562.3730272897010, 66275.9901489112090, 80208.9073262920720, 87936.0963894073648, 122280.763414717993, 130484.415156932562] However, there's still a problem of roundoff error causing the eigenvalues that should be exactly 0 to become nonzero. If (as in this case) the matrix entries are exact rationals, and the matrix is small enough that computing the characteristic polynomial is feasible, you can use a Sturm sequence to avoid trying to compute eigenvalues altogether:

> C:= CharacteristicPolynomial(A, t):
  positives:= sturm(sturmseq(C,t),t,0,infinity);
  zeros:= ldegree(C,t);
  negatives:= degree(C,t)-positives-zeros;
  signature:= positives - negatives;

Actually I think the problem is the "vertical bar" near (but not at) x = 0. This is a consequence of trying to join two points (outside the window) where |cot(x)| is large, one positive and one negative. Try > plot( cot(x), x=-Pi..Pi, view=[-Pi..Pi,-3..3], discont=true );

> PDE := diff(c(x,t),t)=diff(c(x,t),x,x);
  PDE_SOL := pdsolve(PDE);
  ODE1:= op(select(has,op([2,1],PDE_SOL),t));
  ODE2:= op(select(has,op([2,1],PDE_SOL),x));
  dsolve({subs(_c[1]=0,ODE2), _F1(0)=0, _F1(L)=1});
  u := unapply(rhs(%), x);

u := x -> 1/L*x To find the basic solutions with homogeneous boundary conditions:

> _EnvAllSolutions:= true;
  dsolve({ODE2,_F1(0)=0});
  solve(eval(rhs(%),x=L)=0);

{_c[1] = _c[1], _C2 = 0, L = L}, {_c[1] = -Pi^2*_Z1^2/L^2, _C2 = _C2, L = L}

> subs(%[2],%%);

_F1(x) = -_C2*exp((-Pi^2*_Z1^2/L^2)^(1/2)*x)+_C2*exp(-(-Pi^2*_Z1^2/L^2)^(1/2)*x)

> simplify(convert(%,trig)) assuming L>0;

_F1(x) = -2*I*_C2*sin(Pi/L*x*abs(_Z1)) The corresponding solution of the other ODE:

> dsolve(subs(_c[1]=-Pi^2*_Z1^2/L^2,ODE1));

_F2(t) = _C1*exp(-Pi^2*_Z1^2/L^2*t) So the solution of the equation with homogeneous boundary conditions will be

> v := sum(an*sin(Pi*n*x/L)*exp(-Pi^2*n^2*t/L^2),n=1..infinity);

v := sum(an*sin(Pi*n*x/L)*exp(-Pi^2*n^2*t/L^2),n = 1 .. infinity) For t = 0 this is

> eval(v, t=0);

sum(an*sin(Pi*n*x/L),n = 1 .. infinity) This is a Fourier sine series. To find the coefficients:

> an:= 2/L*int((1-u(x))*sin(n*Pi*x/L),x=0..L) assuming n::posint, L > 0;

an := 2/Pi/n So the complete solution is

> u(x) + v;

x/L+sum(2/Pi/n*sin(Pi*n*x/L)*exp(-Pi^2*n^2*t/L^2),n = 1 .. infinity)

Perhaps it should be mentioned that the precise definition of "point of inflection" is not universally agreed upon. Do you include "points" that are not on the curve, i.e. values of x for which f(x) is not defined? Do you include points (c, f(c)) where f'(x) is not defined? And what about cases such as f(x) = x^3 + x^4 sin(1/x), f(0) = 0, where the curve has a tangent line at x=0, is above the tangent line in (0,a) and below it in (-a,0) for some a > 0, but f is not concave up on (0,a) or concave down on (-a,0) for any a > 0? See e.g. the sci.math thread "true or false: (x^5 - x) has inflexion point at origin, (x^4 - x) doesn't" from August 2005, in particular my article of August 24.

Perhaps it should be mentioned that the precise definition of "point of inflection" is not universally agreed upon. Do you include "points" that are not on the curve, i.e. values of x for which f(x) is not defined? Do you include points (c, f(c)) where f'(x) is not defined? And what about cases such as f(x) = x^3 + x^4 sin(1/x), f(0) = 0, where the curve has a tangent line at x=0, is above the tangent line in (0,a) and below it in (-a,0) for some a > 0, but f is not concave up on (0,a) or concave down on (-a,0) for any a > 0? See e.g. the sci.math thread "true or false: (x^5 - x) has inflexion point at origin, (x^4 - x) doesn't" from August 2005, in particular my article of August 24.

The value of the TextArea is a string, namely what is typed in that area.

GetProperty('TextArea0', 'value')

gets that string. Then

parse(...)

takes that string and turns it into a Maple expression (but without evaluating the expression). And finally

SetProperty('MathContainer0', 'value',...)

makes that Maple expression (or rather its typeset version) the value of the Math container. One thing that's annoying: in order to have this take effect, you have to click somewhere in the worksheet or document outside the Text Area and Math Container. You might include a button for people to click on after editing the Text Area. The button doesn't have to actually do anything. If you want to have the expression evaluated instead of just typeset, replace that SetProperty command with

SetProperty('MathContainer0', 'value',
eval(parse(GetProperty('TextArea0', 'value'))));

You could try something like this:

> nestSum:= proc(a,L,r)
   if L = [] then a
   else Sum(nestSum(a,L[2..-1],r),L[1]=r)
   end if
 end;
 Sumpowexpand:= proc(expr)
   local a,v,r,p,L,j; 
   if type(expr,{`*`,`+`}) then 
       return map(procname,expr) 
   end if;
   if not typematch(expr,
            Sum(a::algebraic,v::name=r::range)^p::posint) 
     then return expr 
   end if;
   L:= [seq(v[j],j=1..p)];
   nestSum(mul(subs(v=L[j],a),j=1..p),L,r);
   end;

Then, for example:

> Sumpowexpand(Sum(f(j),j=1..N)^2);

Sum(Sum(f(j[1])*f(j[2]),j[2] = 1 .. N),j[1] = 1 .. N)

You could try something like this:

> nestSum:= proc(a,L,r)
   if L = [] then a
   else Sum(nestSum(a,L[2..-1],r),L[1]=r)
   end if
 end;
 Sumpowexpand:= proc(expr)
   local a,v,r,p,L,j; 
   if type(expr,{`*`,`+`}) then 
       return map(procname,expr) 
   end if;
   if not typematch(expr,
            Sum(a::algebraic,v::name=r::range)^p::posint) 
     then return expr 
   end if;
   L:= [seq(v[j],j=1..p)];
   nestSum(mul(subs(v=L[j],a),j=1..p),L,r);
   end;

Then, for example:

> Sumpowexpand(Sum(f(j),j=1..N)^2);

Sum(Sum(f(j[1])*f(j[2]),j[2] = 1 .. N),j[1] = 1 .. N)

OK, the loop is fine, assuming you have previously assigned values to a, b, i[0], s[0] and r[0]. Now:

If i want to identify the values of S[k], I[k] and R[k] at k=60 how do i enter this?

 
> s[60], i[60], r[60];

And how do i find out the exact point where the epidemic reaches its peak?

I assume that this means the time when i[k] is at its maximum. You could do this:

> imax := max(seq(i[k],k=0..60));
  select(k -> (i[k]=imax), [$0..60]);

OK, the loop is fine, assuming you have previously assigned values to a, b, i[0], s[0] and r[0]. Now:

If i want to identify the values of S[k], I[k] and R[k] at k=60 how do i enter this?

 
> s[60], i[60], r[60];

And how do i find out the exact point where the epidemic reaches its peak?

I assume that this means the time when i[k] is at its maximum. You could do this:

> imax := max(seq(i[k],k=0..60));
  select(k -> (i[k]=imax), [$0..60]);