I would have thought that you could do this by using Worksheet[ReadFile] to get an XML data structure for the worksheet, and then some XMLTools commands to remove the Maple input elements before writing it out again using Worksheet[WriteFile]. This should be just the sort of thing XMLTools is for. But I really haven't used XMLTools, and from looking at the help pages I don't immediately see how to do it. If it is possible, this might be a good thing to include in the help system as an example.

That's not quite how "separation of variables" goes. First of all, the boundary conditions and initial condition are wrong: I'm sure BC2 should be c(L,t) = 1 rather than c(0,L) = 1, and IC should be c(L,0) = 1. Technical point: it's not a good idea to use op(PDE_SOL)[2][1][1] and op(PDE_SOL)[2][1][2] to pick out ODE1 and ODE2: they are in a set, and sets are in arbitrary order, so in a different Maple sessions it's quite possible your ODE1 and ODE2 could be reversed. Better, if you want ODE1 to be the one involving t and ODE2 the one involving x, would be

> ODE1:= op(select(has,op([2,1],PDE_SOL),t));
  ODE2:= op(select(has,op([2,1],PDE_SOL),x));

ODE1 := diff(_F2(t),t) = _c[1]*_F2(t) ODE2 := diff(_F1(x),`$`(x,2)) = _c[1]*_F1(x) 1) Write the solution as c(x,t) = u(x) + v(x,t) where u(x) is a solution of the PDE that is constant in t with boundary conditions u(0) = 0, u(L) = 1. Note that putting in a constant for _F2(t) in ODE1 requires _c[1] = 0. You could use

> dsolve({subs(_c[1]=0,ODE2), _F1(0)=0, _F1(L)=1});
  u := unapply(%, x);

2) Now v(x,t) is a solution of the PDE with homogeneous boundary conditions v(0,t) = v(L,t) = 0 and initial condition v(x,0) = 1 - u(x). Check what nonzero solutions of ODE2 with _F1(0) = 0 satisfy _F1(L) = 0: you'll find _c[1] = -Pi^2*_Z1^2/L^2 (where _Z1 stands for an integer). 3) This corresponds to basic solutions v[n](x,t) = exp(-Pi^2*n^2*t/L^2)*sin(Pi*n*x/L) for positive integers n, and your v(x,t) will be the sum of a series v(x,t) = sum(a[n]*v[n](x,t),n=1..infinity) where the a[n] are chosen to satisfy the initial condition: 1-u(x) = sum(a[n]*sin(Pi*n*x/L), n=1..infinity) That's a Fourier sine series...

I wonder if you're looking at the right function. When a is an odd positive integer, it looks to me like Vq(E,a) -> +infinity as E -> 0+. For example, Vq(E,1) = 1/(12*E) Vq(E,3) = 5/224/E^5-1/(8*E^4)+3/16/E^3 Vq(E,5)=63/5632/E^9-7/64/E^8+105/256/E^7-45/64/E^6+15/32/E^5 On the other hand, when a is not an integer, Vq(E,a) will not be real, because the factor (1-q/2)^(-1+a) will be a non-integer power of a negative number for q > 2.

I wonder if you're looking at the right function. When a is an odd positive integer, it looks to me like Vq(E,a) -> +infinity as E -> 0+. For example, Vq(E,1) = 1/(12*E) Vq(E,3) = 5/224/E^5-1/(8*E^4)+3/16/E^3 Vq(E,5)=63/5632/E^9-7/64/E^8+105/256/E^7-45/64/E^6+15/32/E^5 On the other hand, when a is not an integer, Vq(E,a) will not be real, because the factor (1-q/2)^(-1+a) will be a non-integer power of a negative number for q > 2.

I don't know why you'd want this to return 2 rather than 0. As far as I'm aware, y(t+2) = (something not involving y) should be considered as 0'th order as a difference equation in y. To be of order 2 it would need to involve y(t+2) and y(t) (where t could be the independent variable itself or a translate of it).

I don't know why you'd want this to return 2 rather than 0. As far as I'm aware, y(t+2) = (something not involving y) should be considered as 0'th order as a difference equation in y. To be of order 2 it would need to involve y(t+2) and y(t) (where t could be the independent variable itself or a translate of it).

Yes, that's right. So, to sum up, to have the line y=a*x+b tangent to the curve y=f(x) at x=x1 we need two equations: f(x1)=a*x1+b and f'(x1)=a. To have the same line tangent to the curve y=f(x) at x=x2 we need the same thing but with x1 replaced by x2. So that gives you your four equations in the four unknowns x1, x2, a and b. Hopefully, Maple will be able to solve them. For example:

> f := x -> x^4-6*x^3+13*x^2+x+1;
  equations := {f(x1) = a*x1 + b, f(x2) = a*x2 + b,
    D(f)(x1) = a, D(f)(x2) = a};
  solve(equations);

{x1 = x2, a = 4*x2^3-18*x2^2+26*x2+1, b = -3*x2^4+12*x2^3-13*x2^2+1, x2 = x2}, {b = -3, x2 = 1, a = 13, x1 = 2}, {x1 = 1, b = -3, a = 13, x2 = 2} Ignore the solution with x1 = x2, and take one of the other two (which are the same except for interchanging x1 and x2): x1 = 1 and x2 = 2.

Yes, that's right. So, to sum up, to have the line y=a*x+b tangent to the curve y=f(x) at x=x1 we need two equations: f(x1)=a*x1+b and f'(x1)=a. To have the same line tangent to the curve y=f(x) at x=x2 we need the same thing but with x1 replaced by x2. So that gives you your four equations in the four unknowns x1, x2, a and b. Hopefully, Maple will be able to solve them. For example:

> f := x -> x^4-6*x^3+13*x^2+x+1;
  equations := {f(x1) = a*x1 + b, f(x2) = a*x2 + b,
    D(f)(x1) = a, D(f)(x2) = a};
  solve(equations);

{x1 = x2, a = 4*x2^3-18*x2^2+26*x2+1, b = -3*x2^4+12*x2^3-13*x2^2+1, x2 = x2}, {b = -3, x2 = 1, a = 13, x1 = 2}, {x1 = 1, b = -3, a = 13, x2 = 2} Ignore the solution with x1 = x2, and take one of the other two (which are the same except for interchanging x1 and x2): x1 = 1 and x2 = 2.

"Guaranteed to exist" is not the same thing as "easy to calculate". The continuous antiderivative that the FTC guarantees is the definite integral from some fixed origin to x. But that's not useful if you want a symbolic antiderivative (especially if you'll want to use that symbolic antiderivative to calculate the definite integral). For example, consider f(x) = max(0, g(x)) where g has a known, continuous antiderivative G. So a (discontinuous) antiderivative of f is piecewise(g(x) > 0, G(x), 0). But to construct a continuous antiderivative of f you'd need to know all the points where g changes sign.

For an example of constructing the second axis "by hand", try the following (with DJKeenan's second p12):

> ax2:= plot({[[1,0],[1,100]],
               seq([[1,10*j],[1.015+0.005*(-1)^j,10*j]],
                   j=1..10)},colour=blue),
       seq(textplot([1.02,20*j,0.2*j],colour=blue,align=right),
           j=0..5):
 display([p12, ax2]);

This is not a really satisfactory solution: it requires too much fiddling, and the tickmarks change their lengths when the plot is stretched in the horizontal direction.

For an example of constructing the second axis "by hand", try the following (with DJKeenan's second p12):

> ax2:= plot({[[1,0],[1,100]],
               seq([[1,10*j],[1.015+0.005*(-1)^j,10*j]],
                   j=1..10)},colour=blue),
       seq(textplot([1.02,20*j,0.2*j],colour=blue,align=right),
           j=0..5):
 display([p12, ax2]);

This is not a really satisfactory solution: it requires too much fiddling, and the tickmarks change their lengths when the plot is stretched in the horizontal direction.

We're making progress, I think. You already said that for the line and curve to pass through the same point (at x = x1, let's say) you need f(x1) = a*x1 + b. Now what equation says that they are going "in the same direction", i.e. have the same slope, there?

We're making progress, I think. You already said that for the line and curve to pass through the same point (at x = x1, let's say) you need f(x1) = a*x1 + b. Now what equation says that they are going "in the same direction", i.e. have the same slope, there?

There's no provision in Maple for having two different y axes in the same plot. The closest thing it offers is an array of plots, where you can have your plots side-by-side. Otherwise, you might try constructing what looks like a second axis "by hand" using line segments and textplot.

There's no provision in Maple for having two different y axes in the same plot. The closest thing it offers is an array of plots, where you can have your plots side-by-side. Otherwise, you might try constructing what looks like a second axis "by hand" using line segments and textplot.