Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

This site is on the blink again and won't let me display my worksheet.   Here is the link to a worksheet that shows how to solve the system of PDEs through finite differences.

Download strange-PDE.mw

restart;

with(plots):

A1 := proc(t)
        pointplot([t,t], symbol=solidcircle, symbolsize=50, color=red);
end proc:

A2 := proc(t)
        pointplot([t,t], symbol=solidcircle, symbolsize=50, color=blue);
end proc:

display([
        animate(A1, [t], t=-1..0),
        animate(A2, [t], t=0..1)
], insequence);

 

 

Download mw.mw

 

Your equations are uncoupled.  You may call pdsolve() to solve each of them separately.  That works on all versions of Maple that I have tried.

In principle, pdsolve() should be able to solve the two (uncoupled) equations together as a system, but that fails in Maple 2022 as we see in the attached worksheet.

Here is how things work on solving the system in various versions of Maple:

  • Maple 2017 and 2018 return no solutions, and no errors;

  • Maple 2019, 2020, 2021 return the correct solution;

  • Maple 2022 fails on error.

restart;

kernelopts(version);

`Maple 2022.2, X86 64 LINUX, Oct 23 2022, Build ID 1657361`

Physics:-Version();

`The "Physics Updates" version in the MapleCloud is 1354. The version installed in this computer is 1342 created 2022, November 2, 11:29 hours Pacific Time, found in the directory /home/rouben/maple/toolbox/2022/Physics Updates/lib/`

Pdsolve gets the correct solution here:

pde1 := diff(u(x,t),t) = 0;
ic1 := u(x,0) = f(x);

diff(u(x, t), t) = 0

u(x, 0) = f(x)

pdsolve({pde1,ic1});

u(x, t) = f(x)

Pdsolve gets also the correct solution here:

pde2 := diff(v(x,t),t) = 0;
ic2 := v(x,0) = g(x);

diff(v(x, t), t) = 0

v(x, 0) = g(x)

pdsolve({pde2,ic2});

v(x, t) = g(x)

BUG: Pdsolve gets confused when solving the two (uncoupled) equations together:

pdsolve({pde1, pde2, ic1, ic2});

Error, (in pdsolve) invalid input: indets expects 1 or 2 arguments, but received 3

 

Download mw.mw

I am not going to do the entire homework for you but here is something to get you started.
 

plot(8+8*sin(theta), theta=0..2*Pi, coords=polar, scaling=constrained);

 

plot3d(25 - x^2, x=0..5, y=1..8, view=[0..5,0..8,default]);

 

plots:-display(
    plot3d([10*cos(t)*cos(s), 10*cos(t)*sin(s), 10*sin(t)], t=-arccos(6/10)..arccos(6/10), s=-Pi..Pi),
    plot3d([6*cos(t), 6*sin(t), s], t=-Pi..Pi, s=-8..8),
color=["Green","Red"], style=surface, scaling=constrained, lightmodel=light4);

 


 

Download mw.mw

 

I have not examined the source of the error in your worksheet, but I would suggest using the facilities provided in the Physics[Vectors] package for encoding the Rodrigues formula rather than homogeneous transformation matrices, because the math works more naturally this way.  This worksheet shows how.

restart;

with(Physics[Vectors]):

Rodrigues'  formula:
Rotates the vectorNULL`#mover(mi("v"),mo("→"))` about the unit vecor `#mover(mi("k"),mo("→"))` by the angle alpha:

v__rot_ := v_*cos(alpha) + (k_ &x v_)*sin(alpha) + (k_ . v_)*(1 - cos(alpha))*k_;

v_*cos(alpha)+Physics:-Vectors:-`&x`(k_, v_)*sin(alpha)+Physics:-Vectors:-`.`(k_, v_)*(1-cos(alpha))*k_

Make that into a procedure: Usage:  v__rot = R(`#mover(mi("v"),mo("→"))`, `#mover(mi("k"),mo("→"))`, alpha)

R := unapply(v__rot_, v_, k_, alpha);

proc (v_, k_, alpha) options operator, arrow; v_*cos(alpha)+Physics:-Vectors:-`&x`(k_, v_)*sin(alpha)+Physics:-Vectors:-`.`(k_, v_)*(1-cos(alpha))*k_ end proc

Example 1

 

Pick a unit vector:

1*_i + 2*_j + 3*_k:
k_ := % / sqrt(% . %);

(1/14)*(_i+2*_j+3*_k)*14^(1/2)

Pick an arbitrary vector:

v_ := _i - 3*_j + 2*_k;

_i-3*_j+2*_k

Pick a rotation angle:

alpha := Pi/6;

(1/6)*Pi

Find the rotated vector:

R(v_, k_, alpha):
collect(%, [_i, _j, _k]);

((13/28)*3^(1/2)+(13/28)*14^(1/2)+1/14)*_i+(-(11/7)*3^(1/2)+(1/28)*14^(1/2)+1/7)*_j+((25/28)*3^(1/2)-(5/28)*14^(1/2)+3/14)*_k

 

Example 2

 

Composite rotation

A unit vector:

k__1_ := k_;

(1/14)*(_i+2*_j+3*_k)*14^(1/2)

Another unit vector

k__2_ := subs(_j=-_j, k__1_);

(1/14)*(_i-2*_j+3*_k)*14^(1/2)

Rotate `#mover(mi("v"),mo("→"))` about `#msub(mi("k"),mi("1_"))`by (1/6)*Pi and then rotate the result about `#mi("\`k__2\`")` by (1/4)*Pi:

R(R(v_, k__1_, Pi/6), k__2_, Pi/4):
collect(%, [_i, _j, _k], simplify);

((1/392)*(3*3^(1/2)-4*7^(1/2)+57)*2^(1/2)+(1/392)*(82*7^(1/2)+176)*3^(1/2)+(81/196)*7^(1/2)+3/98)*_i+((1/196)*(-66*3^(1/2)+4*7^(1/2)+174)*2^(1/2)+(1/196)*(7*7^(1/2)-176)*3^(1/2)+(3/196)*7^(1/2)-3/49)*_j+((1/392)*(-89*3^(1/2)-12*7^(1/2)+213)*2^(1/2)+(1/392)*(-18*7^(1/2)+528)*3^(1/2)-(25/196)*7^(1/2)+9/98)*_k

 

 

Download mw.mw

 

 

 

restart;

kernelopts(version);

`Maple 2022.2, X86 64 LINUX, Oct 23 2022, Build ID 1657361`

pd1 := diff(u(x,t),t) = diff(u(x,t),x,x) + (1-alpha)*u(x,t);

diff(u(x, t), t) = diff(diff(u(x, t), x), x)+(1-alpha)*u(x, t)

pd2 := diff(v(x,t),t) = mu*diff(v(x,t),x,x) + beta*v(x,t) + alpha*u(x,t);

diff(v(x, t), t) = mu*(diff(diff(v(x, t), x), x))+beta*v(x, t)+alpha*u(x, t)

ics := u(x,0) = Dirac(x), v(x,0) = 0;

u(x, 0) = Dirac(x), v(x, 0) = 0

Maple is unable solve the system sympolically on its own:

pdsolve({pd1, pd2, ics}) assuming t > 0, mu > 0;

but it can, with a little help.

 

The first PDE has only the unknown u, so we solve it for u:

sol1 := pdsolve({pd1,ics[1]}) assuming t > 0;

u(x, t) = (1/2)*exp((1/4)*((-4*alpha+4)*t^2-x^2)/t)/(t^(1/2)*Pi^(1/2))

Then, substitute the result in the second PDE, and solve that for v:

subs(sol1, pd2);
sol2 := pdsolve({%, ics[2]}) assuming t > 0, mu > 0;

diff(v(x, t), t) = mu*(diff(diff(v(x, t), x), x))+beta*v(x, t)+(1/2)*alpha*exp((1/4)*((-4*alpha+4)*t^2-x^2)/t)/(t^(1/2)*Pi^(1/2))

v(x, t) = alpha*((erf((1/2)*(-2*((-alpha-beta+1)/(mu-1))^(1/2)*t+x)/t^(1/2))-erf((1/2)*(-2*t*mu*((-alpha-beta+1)/(mu-1))^(1/2)+x)/(mu^(1/2)*t^(1/2))))*exp((-x*(mu-1)*((-alpha-beta+1)/(mu-1))^(1/2)-((-1+alpha)*mu+beta)*t)/(mu-1))-exp((x*(mu-1)*((-alpha-beta+1)/(mu-1))^(1/2)-((-1+alpha)*mu+beta)*t)/(mu-1))*(erf((1/2)*(2*((-alpha-beta+1)/(mu-1))^(1/2)*t+x)/t^(1/2))-erf((1/2)*(2*t*mu*((-alpha-beta+1)/(mu-1))^(1/2)+x)/(mu^(1/2)*t^(1/2)))))/(((-alpha-beta+1)/(mu-1))^(1/2)*(4*mu-4))

 

Download mw.mw

 

In several places you have used square brackets [...] and curly braces {...} for grouping your mathematical terms.  In Maple, the only delimiters that you can use for that purpose are the parentheses (...).

Adjust your code as necessary and write back if it still does not work.

delvin, you should know that your equations are not comprehensible to people who are unfamiliar with numerals written in the Farsi/Persian script.  If you desire to reach the widest available help, make an effort to write your mathematics in symbols which most of the readers of this forum can understand.

Anyway, here is the answer to your question.

 

restart;

You wish to perform your hand calculation in Maple.  Here is the original expression:

c := (-exp(4*d)*beta^2 + 4*alpha*exp(4*d) + beta^2 - 4*alpha)/(8*exp(2*d));

(1/8)*(-exp(4*d)*beta^2+4*alpha*exp(4*d)+beta^2-4*alpha)/exp(2*d)

and here is the simplified form obtained by Maple:

convert(simplify(c), trigh);

(-(1/4)*beta^2+alpha)*sinh(2*d)

      That agrees with your hand calculation.

 

 

 

restart;

g := x -> piecewise(x - 1/3 = 0, 1, 0);

g := proc (x) options operator, arrow; piecewise(x-1/3 = 0, 1, 0) end proc

g(3+1/3);

0

g(0+1/3);

1

 

Download mm.mw

 

 

restart;

Take m = 1 and renane "(∂phi)/(∂ xi)" to y:

1/2*y^2 + (1-cos(phi)) = h+2:
solve(%, y):
f := unapply(%[1], phi, h);

proc (phi, h) options operator, arrow; (2*cos(phi)+2*h+2)^(1/2) end proc

seq(f(phi,h), h in [-1,0,1]), seq(-f(phi,h), h in [-1,0,1]);
plot([%], phi=-Pi..3*Pi, size=[800,500]);

2^(1/2)*cos(phi)^(1/2), (2*cos(phi)+2)^(1/2), (2*cos(phi)+4)^(1/2), -2^(1/2)*cos(phi)^(1/2), -(2*cos(phi)+2)^(1/2), -(2*cos(phi)+4)^(1/2)

 

Download mw.mw

 

That initial value problem has infinitely many solutions.  Have a look at this worksheet.

restart;

ode := diff(v(t),t) = -2*v(t)^(2/3);

diff(v(t), t) = -2*v(t)^(2/3)

Find the general solution

dsolve(ode);
dsol := isolate(%, v(t));

v(t)^(1/3)+(2/3)*t-c__1 = 0

v(t) = (-(2/3)*t+c__1)^3

Let c__1 = 0:

eval(dsol, c__1 = 0);

v(t) = -(8/27)*t^3

That's a solution of the ODE with the initial condition v(0) = 0.

That, however, is not the only solution.  The following is also

a solution of the initial value problem for any choice of a > 0NULL

sol := piecewise(t < a, 0, -(8*(t-a)^3)/27);

sol := piecewise(t < a, 0, -8*(t-a)^3*(1/27))

Here is what the solution looks like when a = 2:

eval(sol, a=2);
plot(%, t=0..4, color=red, thickness=5);

piecewise(t < 2, 0, -8*(t-2)^3*(1/27))

Download mw.mw

 

 

The trailing tilde is not a part of  Omega; it is there to remind you that you have made assumptions on it.

If you don't want to see those reminders, disable them with inserting interface(showassumed=0); near the top of your worksheet.

mtaylor does not like expressions with "diff" in them, but it works well when converted to "D".  So all you need is:

convert(lhs(Eq1), D):
mtaylor(%, [x,t], 2);

 

I am posting this as an "Answer" since the question answered here is quite distinct from the question that was originally asked and answered in this thread.

Geodesics connecting two points on a torus

restart;

with(plots):

with(VariationalCalculus):

The parametric equation of a torus with major and minor radii of a and b.

< (a+b*cos(v))*cos(u), (a+b*cos(v))*sin(u),  b*sin(v) >;
x := unapply(%, u, v):

Vector(3, {(1) = (a+b*cos(v))*cos(u), (2) = (a+b*cos(v))*sin(u), (3) = b*sin(v)})

Calculate the coefficients E, F, and G of the first fundamental form:

xu := diff(x(u,v),u):
xv := diff(x(u,v),v):
E := unapply(simplify(xu^+ . xu), u, v);
F := unapply(simplify(xu^+ . xv), u, v);
G := unapply(simplify(xv^+ . xv), u, v);

proc (u, v) options operator, arrow; (a+b*cos(v))^2 end proc

proc (u, v) options operator, arrow; 0 end proc

proc (u, v) options operator, arrow; b^2 end proc

Look for geodesics of the form x(u, v(u)).  These account for all geodesics with the exception

of those corresponding to "u="constant.  The latter are plain circles so we don't need Maple

to calculate them.

 

From differential geometry, the length of  any curve of the form x(u, v(u)), u__1 < u and u < u__2 

is given by L = int(Upsilon(u), u = u__1 .. u__2) where:

Upsilon := sqrt(E(u, v(u)) + 2*F(u, v(u))*diff(v(u),u) + G(u,v(u))*(diff(v(u),u))^2);

((a+b*cos(v(u)))^2+b^2*(diff(v(u), u))^2)^(1/2)

The minimum length is obtained with the help from the Euler-Lagrange equation

which leads to a second order nonlinear differential equation in v(u):

EulerLagrange(Upsilon, u, v(u)):
remove(type, %, equation)[]:
isolate(%, diff(v(u), u, u)):
de := simplify(%);

diff(diff(v(u), u), u) = -(2*b^2*(diff(v(u), u))^2+(a+b*cos(v(u)))^2)*sin(v(u))/((a+b*cos(v(u)))*b)

Pick numbers for the torus's radii:

a, b := 5, 2;

5, 2

torus := plot3d(x(u,v), u=-Pi..Pi, v=-Pi..Pi,
        scaling=constrained, color="khaki", style=surface);

Example 1

 

Calculate a geodesic from u = 0, v = 0 to u = (1/2)*Pi,  v = (1/2)*Pi:

bc1 := v(0)=0,  v(Pi/2)=Pi/2;

v(0) = 0, v((1/2)*Pi) = (1/2)*Pi

dsol1 := dsolve({de, bc1}, numeric, output=operator);

[u = proc (u) local _res, _dat, _solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; _dat := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(21, {(1) = .0, (2) = 0.5860933145427747e-1, (3) = .11721866290855494, (4) = .17582799436283242, (5) = .23445728501980884, (6) = .29400251389637233, (7) = .35502246529707177, (8) = .4176324593904259, (9) = .48228601222277023, (10) = .5512774626633876, (11) = .6260306219686177, (12) = .7081249005896944, (13) = .7985125481567577, (14) = .8998553552891154, (15) = 1.0097294694706644, (16) = 1.1235608288918082, (17) = 1.236518680079433, (18) = 1.3442780572143578, (19) = 1.4354896847639433, (20) = 1.5083320801396727, (21) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(21, 2, {(1, 1) = .0, (1, 2) = 4.357909520225549, (2, 1) = .2533473418001706, (2, 2) = 4.253184219882809, (3, 1) = .4950045392990294, (3, 2) = 3.9682426260835926, (4, 1) = .7163213657367207, (4, 2) = 3.5718304902965725, (5, 1) = .912986520750112, (5, 2) = 3.1352819174892472, (6, 1) = 1.0866622862817927, (6, 2) = 2.7029203986640025, (7, 1) = 1.2390473250060285, (7, 2) = 2.2996023968537105, (8, 1) = 1.3714065523042172, (8, 2) = 1.9372281704148966, (9, 1) = 1.4859944998591417, (9, 2) = 1.6158898922071914, (10, 1) = 1.587137032220874, (10, 2) = 1.3241155689575193, (11, 1) = 1.6758207471908326, (11, 2) = 1.0559025112956482, (12, 1) = 1.7519378227487177, (12, 2) = .8049893321158691, (13, 1) = 1.8136775679360608, (13, 2) = .566571063410118, (14, 1) = 1.8589077269041003, (14, 2) = .33025980575675135, (15, 1) = 1.8821194017143614, (15, 2) = 0.9449920022222437e-1, (16, 1) = 1.879341830166544, (16, 2) = -.1435752401424998, (17, 1) = 1.8494009920056307, (17, 2) = -.3895492919942197, (18, 1) = 1.7937205716265703, (18, 2) = -.649482813775163, (19, 1) = 1.7232026399266667, (19, 2) = -.9032024763979258, (20, 1) = 1.649041988044007, (20, 2) = -1.138744988151322, (21, 1) = 1.5707963267949, (21, 2) = -1.37200094134623}, datatype = float[8], order = C_order); YP := Matrix(21, 2, {(1, 1) = 4.357909520225549, (1, 2) = -.0, (2, 1) = 4.253184219882809, (2, 2) = -3.48400814417279, (3, 1) = 3.9682426260835926, (3, 2) = -6.031898300951247, (4, 1) = 3.5718304902965725, (4, 2) = -7.285198339756213, (5, 1) = 3.1352819174892472, (5, 2) = -7.462358620559608, (6, 1) = 2.7029203986640025, (6, 2) = -6.985675492056729, (7, 1) = 2.2996023968537105, (7, 2) = -6.210452323447571, (8, 1) = 1.9372281704148966, (8, 2) = -5.371378792594195, (9, 1) = 1.6158898922071914, (9, 2) = -4.588579439987948, (10, 1) = 1.3241155689575193, (10, 2) = -3.8949934400656105, (11, 1) = 1.0559025112956482, (11, 2) = -3.3078240571262456, (12, 1) = .8049893321158691, (12, 2) = -2.831415429501529, (13, 1) = .566571063410118, (13, 2) = -2.46897776239132, (14, 1) = .33025980575675135, (14, 2) = -2.2189138501893897, (15, 1) = 0.9449920022222437e-1, (15, 2) = -2.0959796964461206, (16, 1) = -.1435752401424998, (16, 2) = -2.110492861185544, (17, 1) = -.3895492919942197, (17, 2) = -2.2703348506653094, (18, 1) = -.649482813775163, (18, 2) = -2.5835657026971863, (19, 1) = -.9032024763979258, (19, 2) = -3.0077240548867414, (20, 1) = -1.138744988151322, (20, 2) = -3.482022199297429, (21, 1) = -1.37200094134623, (21, 2) = -4.005909266443952}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(21, {(1) = .0, (2) = 0.5860933145427747e-1, (3) = .11721866290855494, (4) = .17582799436283242, (5) = .23445728501980884, (6) = .29400251389637233, (7) = .35502246529707177, (8) = .4176324593904259, (9) = .48228601222277023, (10) = .5512774626633876, (11) = .6260306219686177, (12) = .7081249005896944, (13) = .7985125481567577, (14) = .8998553552891154, (15) = 1.0097294694706644, (16) = 1.1235608288918082, (17) = 1.236518680079433, (18) = 1.3442780572143578, (19) = 1.4354896847639433, (20) = 1.5083320801396727, (21) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(21, 2, {(1, 1) = .0, (1, 2) = 0.6072503241630317e-7, (2, 1) = 0.1633853429871433e-7, (2, 2) = -0.17727261949849427e-7, (3, 1) = 0.61119342794151085e-8, (3, 2) = -0.8326391241592253e-7, (4, 1) = -0.8900614449991116e-8, (4, 2) = -0.5046791702922638e-7, (5, 1) = -0.12191319360321298e-7, (5, 2) = -0.2346020405227589e-7, (6, 1) = -0.9575722148775851e-8, (6, 2) = -0.25118370529491878e-7, (7, 1) = -0.7857479569879596e-8, (7, 2) = -0.28341312692007233e-7, (8, 1) = -0.854262175747852e-8, (8, 2) = -0.22349665722449066e-7, (9, 1) = -0.10412904856270964e-7, (9, 2) = -0.12056403907819441e-7, (10, 1) = -0.12496603030736918e-7, (10, 2) = -0.21735212793440864e-8, (11, 1) = -0.14343296103736601e-7, (11, 2) = 0.4884974297168228e-8, (12, 1) = -0.15809361149223226e-7, (12, 2) = 0.8424061973135961e-8, (13, 1) = -0.1684004558944395e-7, (13, 2) = 0.8468555558744448e-8, (14, 1) = -0.17663088227947315e-7, (14, 2) = 0.5752865457923729e-8, (15, 1) = -0.1799958791726186e-7, (15, 2) = 0.8947994861895889e-9, (16, 1) = -0.1754706929898934e-7, (16, 2) = -0.4029590222151902e-8, (17, 1) = -0.15600567545640277e-7, (17, 2) = -0.5508170053250888e-8, (18, 1) = -0.1101866414537828e-7, (18, 2) = 0.17737187382485602e-8, (19, 1) = -0.6481329226736263e-8, (19, 2) = 0.1388172656806454e-7, (20, 1) = -0.31379335721370494e-8, (20, 2) = 0.2492283863873428e-7, (21, 1) = .0, (21, 2) = 0.3612844847603369e-7}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[21] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(8.326391241592253e-8) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 21, [v(u), diff(v(u), u)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[21] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[21] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(21, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(21, 2, X, Y, outpoint, yout, L, V) end if; [u = outpoint, seq('[v(u), diff(v(u), u)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[21] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(8.326391241592253e-8) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 21, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[21] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[21] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(21, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(21, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 36893628625364698316, (2) = 36893628625364698756, (3) = 36893628625364698932}), (3) = [u, v(u), diff(v(u), u)], (4) = 0}); _solnproc := _dat[1]; if member(u, ["last", 'last']) then _res := _solnproc("last"); if type(_res, 'list') then return _res[1] end if elif type(u, `=`) and member(lhs(u), ["initial", 'initial']) then if type(rhs(u), 'list') then _res := _solnproc("initial" = [0, op(rhs(u))]) else _res := _solnproc("initial" = [1, rhs(u)]) end if; if type(_res, 'list') then return _res[1] end if elif u = "sysvars" then return _dat[3] end if; u end proc, v = proc (u) local res, data, solnproc, v, outpoint; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](u) else outpoint := evalf(u) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(21, {(1) = .0, (2) = 0.5860933145427747e-1, (3) = .11721866290855494, (4) = .17582799436283242, (5) = .23445728501980884, (6) = .29400251389637233, (7) = .35502246529707177, (8) = .4176324593904259, (9) = .48228601222277023, (10) = .5512774626633876, (11) = .6260306219686177, (12) = .7081249005896944, (13) = .7985125481567577, (14) = .8998553552891154, (15) = 1.0097294694706644, (16) = 1.1235608288918082, (17) = 1.236518680079433, (18) = 1.3442780572143578, (19) = 1.4354896847639433, (20) = 1.5083320801396727, (21) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(21, 2, {(1, 1) = .0, (1, 2) = 4.357909520225549, (2, 1) = .2533473418001706, (2, 2) = 4.253184219882809, (3, 1) = .4950045392990294, (3, 2) = 3.9682426260835926, (4, 1) = .7163213657367207, (4, 2) = 3.5718304902965725, (5, 1) = .912986520750112, (5, 2) = 3.1352819174892472, (6, 1) = 1.0866622862817927, (6, 2) = 2.7029203986640025, (7, 1) = 1.2390473250060285, (7, 2) = 2.2996023968537105, (8, 1) = 1.3714065523042172, (8, 2) = 1.9372281704148966, (9, 1) = 1.4859944998591417, (9, 2) = 1.6158898922071914, (10, 1) = 1.587137032220874, (10, 2) = 1.3241155689575193, (11, 1) = 1.6758207471908326, (11, 2) = 1.0559025112956482, (12, 1) = 1.7519378227487177, (12, 2) = .8049893321158691, (13, 1) = 1.8136775679360608, (13, 2) = .566571063410118, (14, 1) = 1.8589077269041003, (14, 2) = .33025980575675135, (15, 1) = 1.8821194017143614, (15, 2) = 0.9449920022222437e-1, (16, 1) = 1.879341830166544, (16, 2) = -.1435752401424998, (17, 1) = 1.8494009920056307, (17, 2) = -.3895492919942197, (18, 1) = 1.7937205716265703, (18, 2) = -.649482813775163, (19, 1) = 1.7232026399266667, (19, 2) = -.9032024763979258, (20, 1) = 1.649041988044007, (20, 2) = -1.138744988151322, (21, 1) = 1.5707963267949, (21, 2) = -1.37200094134623}, datatype = float[8], order = C_order); YP := Matrix(21, 2, {(1, 1) = 4.357909520225549, (1, 2) = -.0, (2, 1) = 4.253184219882809, (2, 2) = -3.48400814417279, (3, 1) = 3.9682426260835926, (3, 2) = -6.031898300951247, (4, 1) = 3.5718304902965725, (4, 2) = -7.285198339756213, (5, 1) = 3.1352819174892472, (5, 2) = -7.462358620559608, (6, 1) = 2.7029203986640025, (6, 2) = -6.985675492056729, (7, 1) = 2.2996023968537105, (7, 2) = -6.210452323447571, (8, 1) = 1.9372281704148966, (8, 2) = -5.371378792594195, (9, 1) = 1.6158898922071914, (9, 2) = -4.588579439987948, (10, 1) = 1.3241155689575193, (10, 2) = -3.8949934400656105, (11, 1) = 1.0559025112956482, (11, 2) = -3.3078240571262456, (12, 1) = .8049893321158691, (12, 2) = -2.831415429501529, (13, 1) = .566571063410118, (13, 2) = -2.46897776239132, (14, 1) = .33025980575675135, (14, 2) = -2.2189138501893897, (15, 1) = 0.9449920022222437e-1, (15, 2) = -2.0959796964461206, (16, 1) = -.1435752401424998, (16, 2) = -2.110492861185544, (17, 1) = -.3895492919942197, (17, 2) = -2.2703348506653094, (18, 1) = -.649482813775163, (18, 2) = -2.5835657026971863, (19, 1) = -.9032024763979258, (19, 2) = -3.0077240548867414, (20, 1) = -1.138744988151322, (20, 2) = -3.482022199297429, (21, 1) = -1.37200094134623, (21, 2) = -4.005909266443952}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(21, {(1) = .0, (2) = 0.5860933145427747e-1, (3) = .11721866290855494, (4) = .17582799436283242, (5) = .23445728501980884, (6) = .29400251389637233, (7) = .35502246529707177, (8) = .4176324593904259, (9) = .48228601222277023, (10) = .5512774626633876, (11) = .6260306219686177, (12) = .7081249005896944, (13) = .7985125481567577, (14) = .8998553552891154, (15) = 1.0097294694706644, (16) = 1.1235608288918082, (17) = 1.236518680079433, (18) = 1.3442780572143578, (19) = 1.4354896847639433, (20) = 1.5083320801396727, (21) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(21, 2, {(1, 1) = .0, (1, 2) = 0.6072503241630317e-7, (2, 1) = 0.1633853429871433e-7, (2, 2) = -0.17727261949849427e-7, (3, 1) = 0.61119342794151085e-8, (3, 2) = -0.8326391241592253e-7, (4, 1) = -0.8900614449991116e-8, (4, 2) = -0.5046791702922638e-7, (5, 1) = -0.12191319360321298e-7, (5, 2) = -0.2346020405227589e-7, (6, 1) = -0.9575722148775851e-8, (6, 2) = -0.25118370529491878e-7, (7, 1) = -0.7857479569879596e-8, (7, 2) = -0.28341312692007233e-7, (8, 1) = -0.854262175747852e-8, (8, 2) = -0.22349665722449066e-7, (9, 1) = -0.10412904856270964e-7, (9, 2) = -0.12056403907819441e-7, (10, 1) = -0.12496603030736918e-7, (10, 2) = -0.21735212793440864e-8, (11, 1) = -0.14343296103736601e-7, (11, 2) = 0.4884974297168228e-8, (12, 1) = -0.15809361149223226e-7, (12, 2) = 0.8424061973135961e-8, (13, 1) = -0.1684004558944395e-7, (13, 2) = 0.8468555558744448e-8, (14, 1) = -0.17663088227947315e-7, (14, 2) = 0.5752865457923729e-8, (15, 1) = -0.1799958791726186e-7, (15, 2) = 0.8947994861895889e-9, (16, 1) = -0.1754706929898934e-7, (16, 2) = -0.4029590222151902e-8, (17, 1) = -0.15600567545640277e-7, (17, 2) = -0.5508170053250888e-8, (18, 1) = -0.1101866414537828e-7, (18, 2) = 0.17737187382485602e-8, (19, 1) = -0.6481329226736263e-8, (19, 2) = 0.1388172656806454e-7, (20, 1) = -0.31379335721370494e-8, (20, 2) = 0.2492283863873428e-7, (21, 1) = .0, (21, 2) = 0.3612844847603369e-7}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[21] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(8.326391241592253e-8) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 21, [v(u), diff(v(u), u)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[21] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[21] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(21, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(21, 2, X, Y, outpoint, yout, L, V) end if; [u = outpoint, seq('[v(u), diff(v(u), u)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[21] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(8.326391241592253e-8) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 21, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[21] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[21] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(21, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(21, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 36893628625364698316, (2) = 36893628625364698756, (3) = 36893628625364698932}), (3) = [u, v(u), diff(v(u), u)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(u) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(u) else v := pointto(data[2][2]); return ('v')(u) end if end if; try res := solnproc(outpoint); res[2] catch: error  end try end proc, D(v) = proc (u) local res, data, solnproc, `D(v)`, outpoint; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](u) else outpoint := evalf(u) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(21, {(1) = .0, (2) = 0.5860933145427747e-1, (3) = .11721866290855494, (4) = .17582799436283242, (5) = .23445728501980884, (6) = .29400251389637233, (7) = .35502246529707177, (8) = .4176324593904259, (9) = .48228601222277023, (10) = .5512774626633876, (11) = .6260306219686177, (12) = .7081249005896944, (13) = .7985125481567577, (14) = .8998553552891154, (15) = 1.0097294694706644, (16) = 1.1235608288918082, (17) = 1.236518680079433, (18) = 1.3442780572143578, (19) = 1.4354896847639433, (20) = 1.5083320801396727, (21) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(21, 2, {(1, 1) = .0, (1, 2) = 4.357909520225549, (2, 1) = .2533473418001706, (2, 2) = 4.253184219882809, (3, 1) = .4950045392990294, (3, 2) = 3.9682426260835926, (4, 1) = .7163213657367207, (4, 2) = 3.5718304902965725, (5, 1) = .912986520750112, (5, 2) = 3.1352819174892472, (6, 1) = 1.0866622862817927, (6, 2) = 2.7029203986640025, (7, 1) = 1.2390473250060285, (7, 2) = 2.2996023968537105, (8, 1) = 1.3714065523042172, (8, 2) = 1.9372281704148966, (9, 1) = 1.4859944998591417, (9, 2) = 1.6158898922071914, (10, 1) = 1.587137032220874, (10, 2) = 1.3241155689575193, (11, 1) = 1.6758207471908326, (11, 2) = 1.0559025112956482, (12, 1) = 1.7519378227487177, (12, 2) = .8049893321158691, (13, 1) = 1.8136775679360608, (13, 2) = .566571063410118, (14, 1) = 1.8589077269041003, (14, 2) = .33025980575675135, (15, 1) = 1.8821194017143614, (15, 2) = 0.9449920022222437e-1, (16, 1) = 1.879341830166544, (16, 2) = -.1435752401424998, (17, 1) = 1.8494009920056307, (17, 2) = -.3895492919942197, (18, 1) = 1.7937205716265703, (18, 2) = -.649482813775163, (19, 1) = 1.7232026399266667, (19, 2) = -.9032024763979258, (20, 1) = 1.649041988044007, (20, 2) = -1.138744988151322, (21, 1) = 1.5707963267949, (21, 2) = -1.37200094134623}, datatype = float[8], order = C_order); YP := Matrix(21, 2, {(1, 1) = 4.357909520225549, (1, 2) = -.0, (2, 1) = 4.253184219882809, (2, 2) = -3.48400814417279, (3, 1) = 3.9682426260835926, (3, 2) = -6.031898300951247, (4, 1) = 3.5718304902965725, (4, 2) = -7.285198339756213, (5, 1) = 3.1352819174892472, (5, 2) = -7.462358620559608, (6, 1) = 2.7029203986640025, (6, 2) = -6.985675492056729, (7, 1) = 2.2996023968537105, (7, 2) = -6.210452323447571, (8, 1) = 1.9372281704148966, (8, 2) = -5.371378792594195, (9, 1) = 1.6158898922071914, (9, 2) = -4.588579439987948, (10, 1) = 1.3241155689575193, (10, 2) = -3.8949934400656105, (11, 1) = 1.0559025112956482, (11, 2) = -3.3078240571262456, (12, 1) = .8049893321158691, (12, 2) = -2.831415429501529, (13, 1) = .566571063410118, (13, 2) = -2.46897776239132, (14, 1) = .33025980575675135, (14, 2) = -2.2189138501893897, (15, 1) = 0.9449920022222437e-1, (15, 2) = -2.0959796964461206, (16, 1) = -.1435752401424998, (16, 2) = -2.110492861185544, (17, 1) = -.3895492919942197, (17, 2) = -2.2703348506653094, (18, 1) = -.649482813775163, (18, 2) = -2.5835657026971863, (19, 1) = -.9032024763979258, (19, 2) = -3.0077240548867414, (20, 1) = -1.138744988151322, (20, 2) = -3.482022199297429, (21, 1) = -1.37200094134623, (21, 2) = -4.005909266443952}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(21, {(1) = .0, (2) = 0.5860933145427747e-1, (3) = .11721866290855494, (4) = .17582799436283242, (5) = .23445728501980884, (6) = .29400251389637233, (7) = .35502246529707177, (8) = .4176324593904259, (9) = .48228601222277023, (10) = .5512774626633876, (11) = .6260306219686177, (12) = .7081249005896944, (13) = .7985125481567577, (14) = .8998553552891154, (15) = 1.0097294694706644, (16) = 1.1235608288918082, (17) = 1.236518680079433, (18) = 1.3442780572143578, (19) = 1.4354896847639433, (20) = 1.5083320801396727, (21) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(21, 2, {(1, 1) = .0, (1, 2) = 0.6072503241630317e-7, (2, 1) = 0.1633853429871433e-7, (2, 2) = -0.17727261949849427e-7, (3, 1) = 0.61119342794151085e-8, (3, 2) = -0.8326391241592253e-7, (4, 1) = -0.8900614449991116e-8, (4, 2) = -0.5046791702922638e-7, (5, 1) = -0.12191319360321298e-7, (5, 2) = -0.2346020405227589e-7, (6, 1) = -0.9575722148775851e-8, (6, 2) = -0.25118370529491878e-7, (7, 1) = -0.7857479569879596e-8, (7, 2) = -0.28341312692007233e-7, (8, 1) = -0.854262175747852e-8, (8, 2) = -0.22349665722449066e-7, (9, 1) = -0.10412904856270964e-7, (9, 2) = -0.12056403907819441e-7, (10, 1) = -0.12496603030736918e-7, (10, 2) = -0.21735212793440864e-8, (11, 1) = -0.14343296103736601e-7, (11, 2) = 0.4884974297168228e-8, (12, 1) = -0.15809361149223226e-7, (12, 2) = 0.8424061973135961e-8, (13, 1) = -0.1684004558944395e-7, (13, 2) = 0.8468555558744448e-8, (14, 1) = -0.17663088227947315e-7, (14, 2) = 0.5752865457923729e-8, (15, 1) = -0.1799958791726186e-7, (15, 2) = 0.8947994861895889e-9, (16, 1) = -0.1754706929898934e-7, (16, 2) = -0.4029590222151902e-8, (17, 1) = -0.15600567545640277e-7, (17, 2) = -0.5508170053250888e-8, (18, 1) = -0.1101866414537828e-7, (18, 2) = 0.17737187382485602e-8, (19, 1) = -0.6481329226736263e-8, (19, 2) = 0.1388172656806454e-7, (20, 1) = -0.31379335721370494e-8, (20, 2) = 0.2492283863873428e-7, (21, 1) = .0, (21, 2) = 0.3612844847603369e-7}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[21] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(8.326391241592253e-8) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 21, [v(u), diff(v(u), u)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[21] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[21] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(21, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(21, 2, X, Y, outpoint, yout, L, V) end if; [u = outpoint, seq('[v(u), diff(v(u), u)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[21] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(8.326391241592253e-8) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 21, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[21] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[21] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(21, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(21, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 36893628625364698316, (2) = 36893628625364698756, (3) = 36893628625364698932}), (3) = [u, v(u), diff(v(u), u)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(u) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(u) else `D(v)` := pointto(data[2][3]); return ('`D(v)`')(u) end if end if; try res := solnproc(outpoint); res[3] catch: error  end try end proc]

display(
        torus,
        tubeplot([seq(eval(x(u, v(u)), dsol1))], u=0..Pi/2, color=red, radius=0.1),
        pointplot3d(eval([x(0,v(0)),x(Pi/2,v(Pi/2))], dsol1), symbol=solidsphere, symbolsize=20, color="Orange"),
lightmodel=light4, thickness=5, orientation=[-5,75,0], size=[700,400],
style=surface, axes=none);

 

Example 2

 

bc2 := v(0)=0,  v(Pi/2)=Pi;

v(0) = 0, v((1/2)*Pi) = Pi

dsol2 := dsolve({de, bc2}, numeric, output=operator);

[u = proc (u) local _res, _dat, _solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; _dat := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(25, {(1) = .0, (2) = 0.3556697983665059e-1, (3) = 0.7150815650209673e-1, (4) = .1078355952208751, (5) = .14456202386531472, (6) = .18170088505495013, (7) = .2192663936438809, (8) = .25772582968519775, (9) = .2979698369946872, (10) = .34028667559059855, (11) = .38503401227598166, (12) = .43268598844062967, (13) = .48418267682069066, (14) = .54178398151824, (15) = .6084908129577877, (16) = .6905705623174822, (17) = .7936582413926776, (18) = .9362490342852545, (19) = 1.1290955386275305, (20) = 1.2785289405368818, (21) = 1.3762683057594027, (22) = 1.4410821559755675, (23) = 1.4911439441371346, (24) = 1.5334604576293163, (25) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(25, 2, {(1, 1) = .0, (1, 2) = 8.027222261980723, (2, 1) = .2831190340450847, (2, 2) = 7.828369059489753, (3, 1) = .5555864658560054, (3, 2) = 7.287670355233425, (4, 1) = .8071856825371913, (4, 2) = 6.54481306582631, (5, 1) = 1.0328331222442264, (5, 2) = 5.745156668353998, (6, 1) = 1.2318715807045222, (6, 2) = 4.987033344028576, (7, 1) = 1.4062899344379456, (7, 2) = 4.316682182254274, (8, 1) = 1.5608590578898303, (8, 2) = 3.739193866869814, (9, 1) = 1.7009934741789976, (9, 2) = 3.241862801797417, (10, 1) = 1.8288761198346395, (10, 2) = 2.817246508505084, (11, 1) = 1.9465497589776795, (11, 2) = 2.455373052357135, (12, 1) = 2.0559118912274252, (12, 2) = 2.146059224600939, (13, 1) = 2.1592833618063074, (13, 2) = 1.8786440703412244, (14, 1) = 2.260381800745741, (14, 2) = 1.6409272186274388, (15, 1) = 2.362338845054942, (15, 2) = 1.4249780044463283, (16, 1) = 2.4705572534383973, (16, 2) = 1.2215356731044336, (17, 1) = 2.5862348892169664, (17, 2) = 1.0329309045427477, (18, 1) = 2.7197502962411266, (18, 2) = .8522354096629391, (19, 1) = 2.868097308964122, (19, 2) = .7002166019075773, (20, 1) = 2.9670989492898685, (20, 2) = .6306046195751053, (21, 1) = 3.0272505993270076, (21, 2) = .6023260373483883, (22, 1) = 3.065868563957721, (22, 2) = .5901748294236376, (23, 1) = 3.095251805141905, (23, 2) = .5841893523046803, (24, 1) = 3.11990600006039, (24, 2) = .5813793052108006, (25, 1) = 3.14159265358979, (25, 2) = .5805900104831727}, datatype = float[8], order = C_order); YP := Matrix(25, 2, {(1, 1) = 8.027222261980723, (1, 2) = -.0, (2, 1) = 7.828369059489753, (2, 2) = -10.861805181334308, (3, 1) = 7.287670355233425, (3, 2) = -18.492634875073744, (4, 1) = 6.54481306582631, (4, 2) = -21.694990694475745, (5, 1) = 5.745156668353998, (5, 2) = -21.40568916672294, (6, 1) = 4.987033344028576, (6, 2) = -19.233324103734176, (7, 1) = 4.316682182254274, (7, 2) = -16.42945105476428, (8, 1) = 3.739193866869814, (8, 2) = -13.650236112103359, (9, 1) = 3.241862801797417, (9, 2) = -11.143337685208598, (10, 1) = 2.817246508505084, (10, 2) = -9.007662429827185, (11, 1) = 2.455373052357135, (11, 2) = -7.242681944771625, (12, 1) = 2.146059224600939, (12, 2) = -5.805737354677115, (13, 1) = 1.8786440703412244, (13, 2) = -4.636510741573245, (14, 1) = 1.6409272186274388, (14, 2) = -3.6671563351547207, (15, 1) = 1.4249780044463283, (15, 2) = -2.852577215638926, (16, 1) = 1.2215356731044336, (16, 2) = -2.1483661147204827, (17, 1) = 1.0329309045427477, (17, 2) = -1.5518631600869446, (18, 1) = .8522354096629391, (18, 2) = -1.024669368486231, (19, 1) = .7002166019075773, (19, 2) = -.5874907987071298, (20, 1) = .6306046195751053, (20, 2) = -.35417872498610686, (21, 1) = .6023260373483883, (21, 2) = -.22683550784375403, (22, 1) = .5901748294236376, (22, 2) = -.14876072185358538, (23, 1) = .5841893523046803, (23, 2) = -0.9060027226914852e-1, (24, 1) = .5813793052108006, (24, 2) = -0.4230373846316707e-1, (25, 1) = .5805900104831727, (25, 2) = -0.6298834938184995e-14}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(25, {(1) = .0, (2) = 0.3556697983665059e-1, (3) = 0.7150815650209673e-1, (4) = .1078355952208751, (5) = .14456202386531472, (6) = .18170088505495013, (7) = .2192663936438809, (8) = .25772582968519775, (9) = .2979698369946872, (10) = .34028667559059855, (11) = .38503401227598166, (12) = .43268598844062967, (13) = .48418267682069066, (14) = .54178398151824, (15) = .6084908129577877, (16) = .6905705623174822, (17) = .7936582413926776, (18) = .9362490342852545, (19) = 1.1290955386275305, (20) = 1.2785289405368818, (21) = 1.3762683057594027, (22) = 1.4410821559755675, (23) = 1.4911439441371346, (24) = 1.5334604576293163, (25) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(25, 2, {(1, 1) = .0, (1, 2) = 0.31500831778913094e-6, (2, 1) = 0.3132220020535465e-7, (2, 2) = 0.7779131439750735e-7, (3, 1) = 0.14325613236197187e-7, (3, 2) = -0.6898540876289522e-7, (4, 1) = -0.6036355453398498e-8, (4, 2) = 0.7943001698758059e-7, (5, 1) = -0.13260007889026392e-8, (5, 2) = 0.11756415179265621e-6, (6, 1) = 0.1071471554168767e-7, (6, 2) = 0.44460301809323215e-7, (7, 1) = 0.16948724685756332e-7, (7, 2) = -0.6994714895790386e-8, (8, 1) = 0.1782055631368677e-7, (8, 2) = -0.1641681279404671e-7, (9, 1) = 0.16393950975769278e-7, (9, 2) = -0.58427989726364065e-8, (10, 1) = 0.14618894539761003e-7, (10, 2) = 0.7611526601257115e-8, (11, 1) = 0.13363200796132462e-7, (11, 2) = 0.16937803808645205e-7, (12, 1) = 0.12801606379592921e-7, (12, 2) = 0.21469959634400602e-7, (13, 1) = 0.12794426509551367e-7, (13, 2) = 0.22808086897112244e-7, (14, 1) = 0.13080732166085607e-7, (14, 2) = 0.22688679611190916e-7, (15, 1) = 0.13483594578904676e-7, (15, 2) = 0.2185857499871233e-7, (16, 1) = 0.13492144691744485e-7, (16, 2) = 0.2072575098426527e-7, (17, 1) = 0.12840479015099201e-7, (17, 2) = 0.17970028302554435e-7, (18, 1) = 0.720289485656478e-8, (18, 2) = 0.11283105136137962e-7, (19, 1) = -0.10626354523447358e-8, (19, 2) = -0.7464659107819674e-8, (20, 1) = 0.19375865427224214e-8, (20, 2) = -0.1047570220567349e-7, (21, 1) = 0.15016966694503165e-8, (21, 2) = -0.914149161086985e-8, (22, 1) = 0.10081859027200914e-8, (22, 2) = -0.848749960193092e-8, (23, 1) = 0.618147262744349e-9, (23, 2) = -0.8172877515951791e-8, (24, 1) = 0.2897545785236149e-9, (24, 2) = -0.8027811845254448e-8, (25, 1) = .0, (25, 2) = -0.7987456326182811e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[25] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.1500831778913094e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 25, [v(u), diff(v(u), u)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[25] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[25] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(25, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(25, 2, X, Y, outpoint, yout, L, V) end if; [u = outpoint, seq('[v(u), diff(v(u), u)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[25] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.1500831778913094e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 25, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[25] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[25] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(25, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(25, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 36893628625364672284, (2) = 36893628625364664316, (3) = 36893628625364664492}), (3) = [u, v(u), diff(v(u), u)], (4) = 0}); _solnproc := _dat[1]; if member(u, ["last", 'last']) then _res := _solnproc("last"); if type(_res, 'list') then return _res[1] end if elif type(u, `=`) and member(lhs(u), ["initial", 'initial']) then if type(rhs(u), 'list') then _res := _solnproc("initial" = [0, op(rhs(u))]) else _res := _solnproc("initial" = [1, rhs(u)]) end if; if type(_res, 'list') then return _res[1] end if elif u = "sysvars" then return _dat[3] end if; u end proc, v = proc (u) local res, data, solnproc, v, outpoint; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](u) else outpoint := evalf(u) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(25, {(1) = .0, (2) = 0.3556697983665059e-1, (3) = 0.7150815650209673e-1, (4) = .1078355952208751, (5) = .14456202386531472, (6) = .18170088505495013, (7) = .2192663936438809, (8) = .25772582968519775, (9) = .2979698369946872, (10) = .34028667559059855, (11) = .38503401227598166, (12) = .43268598844062967, (13) = .48418267682069066, (14) = .54178398151824, (15) = .6084908129577877, (16) = .6905705623174822, (17) = .7936582413926776, (18) = .9362490342852545, (19) = 1.1290955386275305, (20) = 1.2785289405368818, (21) = 1.3762683057594027, (22) = 1.4410821559755675, (23) = 1.4911439441371346, (24) = 1.5334604576293163, (25) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(25, 2, {(1, 1) = .0, (1, 2) = 8.027222261980723, (2, 1) = .2831190340450847, (2, 2) = 7.828369059489753, (3, 1) = .5555864658560054, (3, 2) = 7.287670355233425, (4, 1) = .8071856825371913, (4, 2) = 6.54481306582631, (5, 1) = 1.0328331222442264, (5, 2) = 5.745156668353998, (6, 1) = 1.2318715807045222, (6, 2) = 4.987033344028576, (7, 1) = 1.4062899344379456, (7, 2) = 4.316682182254274, (8, 1) = 1.5608590578898303, (8, 2) = 3.739193866869814, (9, 1) = 1.7009934741789976, (9, 2) = 3.241862801797417, (10, 1) = 1.8288761198346395, (10, 2) = 2.817246508505084, (11, 1) = 1.9465497589776795, (11, 2) = 2.455373052357135, (12, 1) = 2.0559118912274252, (12, 2) = 2.146059224600939, (13, 1) = 2.1592833618063074, (13, 2) = 1.8786440703412244, (14, 1) = 2.260381800745741, (14, 2) = 1.6409272186274388, (15, 1) = 2.362338845054942, (15, 2) = 1.4249780044463283, (16, 1) = 2.4705572534383973, (16, 2) = 1.2215356731044336, (17, 1) = 2.5862348892169664, (17, 2) = 1.0329309045427477, (18, 1) = 2.7197502962411266, (18, 2) = .8522354096629391, (19, 1) = 2.868097308964122, (19, 2) = .7002166019075773, (20, 1) = 2.9670989492898685, (20, 2) = .6306046195751053, (21, 1) = 3.0272505993270076, (21, 2) = .6023260373483883, (22, 1) = 3.065868563957721, (22, 2) = .5901748294236376, (23, 1) = 3.095251805141905, (23, 2) = .5841893523046803, (24, 1) = 3.11990600006039, (24, 2) = .5813793052108006, (25, 1) = 3.14159265358979, (25, 2) = .5805900104831727}, datatype = float[8], order = C_order); YP := Matrix(25, 2, {(1, 1) = 8.027222261980723, (1, 2) = -.0, (2, 1) = 7.828369059489753, (2, 2) = -10.861805181334308, (3, 1) = 7.287670355233425, (3, 2) = -18.492634875073744, (4, 1) = 6.54481306582631, (4, 2) = -21.694990694475745, (5, 1) = 5.745156668353998, (5, 2) = -21.40568916672294, (6, 1) = 4.987033344028576, (6, 2) = -19.233324103734176, (7, 1) = 4.316682182254274, (7, 2) = -16.42945105476428, (8, 1) = 3.739193866869814, (8, 2) = -13.650236112103359, (9, 1) = 3.241862801797417, (9, 2) = -11.143337685208598, (10, 1) = 2.817246508505084, (10, 2) = -9.007662429827185, (11, 1) = 2.455373052357135, (11, 2) = -7.242681944771625, (12, 1) = 2.146059224600939, (12, 2) = -5.805737354677115, (13, 1) = 1.8786440703412244, (13, 2) = -4.636510741573245, (14, 1) = 1.6409272186274388, (14, 2) = -3.6671563351547207, (15, 1) = 1.4249780044463283, (15, 2) = -2.852577215638926, (16, 1) = 1.2215356731044336, (16, 2) = -2.1483661147204827, (17, 1) = 1.0329309045427477, (17, 2) = -1.5518631600869446, (18, 1) = .8522354096629391, (18, 2) = -1.024669368486231, (19, 1) = .7002166019075773, (19, 2) = -.5874907987071298, (20, 1) = .6306046195751053, (20, 2) = -.35417872498610686, (21, 1) = .6023260373483883, (21, 2) = -.22683550784375403, (22, 1) = .5901748294236376, (22, 2) = -.14876072185358538, (23, 1) = .5841893523046803, (23, 2) = -0.9060027226914852e-1, (24, 1) = .5813793052108006, (24, 2) = -0.4230373846316707e-1, (25, 1) = .5805900104831727, (25, 2) = -0.6298834938184995e-14}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(25, {(1) = .0, (2) = 0.3556697983665059e-1, (3) = 0.7150815650209673e-1, (4) = .1078355952208751, (5) = .14456202386531472, (6) = .18170088505495013, (7) = .2192663936438809, (8) = .25772582968519775, (9) = .2979698369946872, (10) = .34028667559059855, (11) = .38503401227598166, (12) = .43268598844062967, (13) = .48418267682069066, (14) = .54178398151824, (15) = .6084908129577877, (16) = .6905705623174822, (17) = .7936582413926776, (18) = .9362490342852545, (19) = 1.1290955386275305, (20) = 1.2785289405368818, (21) = 1.3762683057594027, (22) = 1.4410821559755675, (23) = 1.4911439441371346, (24) = 1.5334604576293163, (25) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(25, 2, {(1, 1) = .0, (1, 2) = 0.31500831778913094e-6, (2, 1) = 0.3132220020535465e-7, (2, 2) = 0.7779131439750735e-7, (3, 1) = 0.14325613236197187e-7, (3, 2) = -0.6898540876289522e-7, (4, 1) = -0.6036355453398498e-8, (4, 2) = 0.7943001698758059e-7, (5, 1) = -0.13260007889026392e-8, (5, 2) = 0.11756415179265621e-6, (6, 1) = 0.1071471554168767e-7, (6, 2) = 0.44460301809323215e-7, (7, 1) = 0.16948724685756332e-7, (7, 2) = -0.6994714895790386e-8, (8, 1) = 0.1782055631368677e-7, (8, 2) = -0.1641681279404671e-7, (9, 1) = 0.16393950975769278e-7, (9, 2) = -0.58427989726364065e-8, (10, 1) = 0.14618894539761003e-7, (10, 2) = 0.7611526601257115e-8, (11, 1) = 0.13363200796132462e-7, (11, 2) = 0.16937803808645205e-7, (12, 1) = 0.12801606379592921e-7, (12, 2) = 0.21469959634400602e-7, (13, 1) = 0.12794426509551367e-7, (13, 2) = 0.22808086897112244e-7, (14, 1) = 0.13080732166085607e-7, (14, 2) = 0.22688679611190916e-7, (15, 1) = 0.13483594578904676e-7, (15, 2) = 0.2185857499871233e-7, (16, 1) = 0.13492144691744485e-7, (16, 2) = 0.2072575098426527e-7, (17, 1) = 0.12840479015099201e-7, (17, 2) = 0.17970028302554435e-7, (18, 1) = 0.720289485656478e-8, (18, 2) = 0.11283105136137962e-7, (19, 1) = -0.10626354523447358e-8, (19, 2) = -0.7464659107819674e-8, (20, 1) = 0.19375865427224214e-8, (20, 2) = -0.1047570220567349e-7, (21, 1) = 0.15016966694503165e-8, (21, 2) = -0.914149161086985e-8, (22, 1) = 0.10081859027200914e-8, (22, 2) = -0.848749960193092e-8, (23, 1) = 0.618147262744349e-9, (23, 2) = -0.8172877515951791e-8, (24, 1) = 0.2897545785236149e-9, (24, 2) = -0.8027811845254448e-8, (25, 1) = .0, (25, 2) = -0.7987456326182811e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[25] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.1500831778913094e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 25, [v(u), diff(v(u), u)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[25] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[25] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(25, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(25, 2, X, Y, outpoint, yout, L, V) end if; [u = outpoint, seq('[v(u), diff(v(u), u)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[25] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.1500831778913094e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 25, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[25] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[25] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(25, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(25, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 36893628625364672284, (2) = 36893628625364664316, (3) = 36893628625364664492}), (3) = [u, v(u), diff(v(u), u)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(u) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(u) else v := pointto(data[2][2]); return ('v')(u) end if end if; try res := solnproc(outpoint); res[2] catch: error  end try end proc, D(v) = proc (u) local res, data, solnproc, `D(v)`, outpoint; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](u) else outpoint := evalf(u) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(25, {(1) = .0, (2) = 0.3556697983665059e-1, (3) = 0.7150815650209673e-1, (4) = .1078355952208751, (5) = .14456202386531472, (6) = .18170088505495013, (7) = .2192663936438809, (8) = .25772582968519775, (9) = .2979698369946872, (10) = .34028667559059855, (11) = .38503401227598166, (12) = .43268598844062967, (13) = .48418267682069066, (14) = .54178398151824, (15) = .6084908129577877, (16) = .6905705623174822, (17) = .7936582413926776, (18) = .9362490342852545, (19) = 1.1290955386275305, (20) = 1.2785289405368818, (21) = 1.3762683057594027, (22) = 1.4410821559755675, (23) = 1.4911439441371346, (24) = 1.5334604576293163, (25) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(25, 2, {(1, 1) = .0, (1, 2) = 8.027222261980723, (2, 1) = .2831190340450847, (2, 2) = 7.828369059489753, (3, 1) = .5555864658560054, (3, 2) = 7.287670355233425, (4, 1) = .8071856825371913, (4, 2) = 6.54481306582631, (5, 1) = 1.0328331222442264, (5, 2) = 5.745156668353998, (6, 1) = 1.2318715807045222, (6, 2) = 4.987033344028576, (7, 1) = 1.4062899344379456, (7, 2) = 4.316682182254274, (8, 1) = 1.5608590578898303, (8, 2) = 3.739193866869814, (9, 1) = 1.7009934741789976, (9, 2) = 3.241862801797417, (10, 1) = 1.8288761198346395, (10, 2) = 2.817246508505084, (11, 1) = 1.9465497589776795, (11, 2) = 2.455373052357135, (12, 1) = 2.0559118912274252, (12, 2) = 2.146059224600939, (13, 1) = 2.1592833618063074, (13, 2) = 1.8786440703412244, (14, 1) = 2.260381800745741, (14, 2) = 1.6409272186274388, (15, 1) = 2.362338845054942, (15, 2) = 1.4249780044463283, (16, 1) = 2.4705572534383973, (16, 2) = 1.2215356731044336, (17, 1) = 2.5862348892169664, (17, 2) = 1.0329309045427477, (18, 1) = 2.7197502962411266, (18, 2) = .8522354096629391, (19, 1) = 2.868097308964122, (19, 2) = .7002166019075773, (20, 1) = 2.9670989492898685, (20, 2) = .6306046195751053, (21, 1) = 3.0272505993270076, (21, 2) = .6023260373483883, (22, 1) = 3.065868563957721, (22, 2) = .5901748294236376, (23, 1) = 3.095251805141905, (23, 2) = .5841893523046803, (24, 1) = 3.11990600006039, (24, 2) = .5813793052108006, (25, 1) = 3.14159265358979, (25, 2) = .5805900104831727}, datatype = float[8], order = C_order); YP := Matrix(25, 2, {(1, 1) = 8.027222261980723, (1, 2) = -.0, (2, 1) = 7.828369059489753, (2, 2) = -10.861805181334308, (3, 1) = 7.287670355233425, (3, 2) = -18.492634875073744, (4, 1) = 6.54481306582631, (4, 2) = -21.694990694475745, (5, 1) = 5.745156668353998, (5, 2) = -21.40568916672294, (6, 1) = 4.987033344028576, (6, 2) = -19.233324103734176, (7, 1) = 4.316682182254274, (7, 2) = -16.42945105476428, (8, 1) = 3.739193866869814, (8, 2) = -13.650236112103359, (9, 1) = 3.241862801797417, (9, 2) = -11.143337685208598, (10, 1) = 2.817246508505084, (10, 2) = -9.007662429827185, (11, 1) = 2.455373052357135, (11, 2) = -7.242681944771625, (12, 1) = 2.146059224600939, (12, 2) = -5.805737354677115, (13, 1) = 1.8786440703412244, (13, 2) = -4.636510741573245, (14, 1) = 1.6409272186274388, (14, 2) = -3.6671563351547207, (15, 1) = 1.4249780044463283, (15, 2) = -2.852577215638926, (16, 1) = 1.2215356731044336, (16, 2) = -2.1483661147204827, (17, 1) = 1.0329309045427477, (17, 2) = -1.5518631600869446, (18, 1) = .8522354096629391, (18, 2) = -1.024669368486231, (19, 1) = .7002166019075773, (19, 2) = -.5874907987071298, (20, 1) = .6306046195751053, (20, 2) = -.35417872498610686, (21, 1) = .6023260373483883, (21, 2) = -.22683550784375403, (22, 1) = .5901748294236376, (22, 2) = -.14876072185358538, (23, 1) = .5841893523046803, (23, 2) = -0.9060027226914852e-1, (24, 1) = .5813793052108006, (24, 2) = -0.4230373846316707e-1, (25, 1) = .5805900104831727, (25, 2) = -0.6298834938184995e-14}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(25, {(1) = .0, (2) = 0.3556697983665059e-1, (3) = 0.7150815650209673e-1, (4) = .1078355952208751, (5) = .14456202386531472, (6) = .18170088505495013, (7) = .2192663936438809, (8) = .25772582968519775, (9) = .2979698369946872, (10) = .34028667559059855, (11) = .38503401227598166, (12) = .43268598844062967, (13) = .48418267682069066, (14) = .54178398151824, (15) = .6084908129577877, (16) = .6905705623174822, (17) = .7936582413926776, (18) = .9362490342852545, (19) = 1.1290955386275305, (20) = 1.2785289405368818, (21) = 1.3762683057594027, (22) = 1.4410821559755675, (23) = 1.4911439441371346, (24) = 1.5334604576293163, (25) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(25, 2, {(1, 1) = .0, (1, 2) = 0.31500831778913094e-6, (2, 1) = 0.3132220020535465e-7, (2, 2) = 0.7779131439750735e-7, (3, 1) = 0.14325613236197187e-7, (3, 2) = -0.6898540876289522e-7, (4, 1) = -0.6036355453398498e-8, (4, 2) = 0.7943001698758059e-7, (5, 1) = -0.13260007889026392e-8, (5, 2) = 0.11756415179265621e-6, (6, 1) = 0.1071471554168767e-7, (6, 2) = 0.44460301809323215e-7, (7, 1) = 0.16948724685756332e-7, (7, 2) = -0.6994714895790386e-8, (8, 1) = 0.1782055631368677e-7, (8, 2) = -0.1641681279404671e-7, (9, 1) = 0.16393950975769278e-7, (9, 2) = -0.58427989726364065e-8, (10, 1) = 0.14618894539761003e-7, (10, 2) = 0.7611526601257115e-8, (11, 1) = 0.13363200796132462e-7, (11, 2) = 0.16937803808645205e-7, (12, 1) = 0.12801606379592921e-7, (12, 2) = 0.21469959634400602e-7, (13, 1) = 0.12794426509551367e-7, (13, 2) = 0.22808086897112244e-7, (14, 1) = 0.13080732166085607e-7, (14, 2) = 0.22688679611190916e-7, (15, 1) = 0.13483594578904676e-7, (15, 2) = 0.2185857499871233e-7, (16, 1) = 0.13492144691744485e-7, (16, 2) = 0.2072575098426527e-7, (17, 1) = 0.12840479015099201e-7, (17, 2) = 0.17970028302554435e-7, (18, 1) = 0.720289485656478e-8, (18, 2) = 0.11283105136137962e-7, (19, 1) = -0.10626354523447358e-8, (19, 2) = -0.7464659107819674e-8, (20, 1) = 0.19375865427224214e-8, (20, 2) = -0.1047570220567349e-7, (21, 1) = 0.15016966694503165e-8, (21, 2) = -0.914149161086985e-8, (22, 1) = 0.10081859027200914e-8, (22, 2) = -0.848749960193092e-8, (23, 1) = 0.618147262744349e-9, (23, 2) = -0.8172877515951791e-8, (24, 1) = 0.2897545785236149e-9, (24, 2) = -0.8027811845254448e-8, (25, 1) = .0, (25, 2) = -0.7987456326182811e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[25] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.1500831778913094e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 25, [v(u), diff(v(u), u)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[25] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[25] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(25, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(25, 2, X, Y, outpoint, yout, L, V) end if; [u = outpoint, seq('[v(u), diff(v(u), u)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[25] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(3.1500831778913094e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 25, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[25] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[25] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(25, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(25, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 36893628625364672284, (2) = 36893628625364664316, (3) = 36893628625364664492}), (3) = [u, v(u), diff(v(u), u)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(u) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(u) else `D(v)` := pointto(data[2][3]); return ('`D(v)`')(u) end if end if; try res := solnproc(outpoint); res[3] catch: error  end try end proc]

display(
        torus,
        tubeplot([seq(eval(x(u, v(u)), dsol2))], u=0..Pi/2, color=red, radius=0.1),
        pointplot3d(eval([x(0,v(0)),x(Pi/2,v(Pi/2))], dsol2), symbol=solidsphere, symbolsize=20, color="Orange"),
lightmodel=light4, thickness=5, orientation=[-65,65,0], size=[700,400],
style=surface, axes=none);

 

 

Example 3

 

The end points in this example are identical to those of Example 1 but the new geodesic,
drawn in cyan, is different.  In fact, there are infinitely many  geodesics that connect

those two points.

bc3 := v(0)=0,  v(Pi/2)=5*Pi/2;

v(0) = 0, v((1/2)*Pi) = (5/2)*Pi

dsol3 := dsolve({de, bc3}, numeric, output=operator);

[u = proc (u) local _res, _dat, _solnproc; option `Copyright (c) 1993 by the University of Waterloo. All rights reserved.`; _dat := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(59, {(1) = .0, (2) = 0.13680812436809657e-1, (3) = 0.27971051145644167e-1, (4) = 0.4296023135974479e-1, (5) = 0.5851786182852681e-1, (6) = 0.7435973609701907e-1, (7) = 0.9049963296106847e-1, (8) = .1070194360291177, (9) = .12400404122395811, (10) = .14149911102616744, (11) = .15962080170244286, (12) = .17846661163773977, (13) = .1981659899247577, (14) = .2189484735830971, (15) = .2410342106109193, (16) = .2648591917099649, (17) = .29096256417430655, (18) = .3201405168424954, (19) = .3537506263213988, (20) = .39359925436482457, (21) = .44328213017143814, (22) = .508829664760916, (23) = .5993459008179163, (24) = .7156934169226271, (25) = .8282966673508598, (26) = .9130409367878883, (27) = .9748398788122531, (28) = 1.0223008015052495, (29) = 1.0606661846616776, (30) = 1.0932188074676115, (31) = 1.1216707736458404, (32) = 1.1472399927159063, (33) = 1.1706725764522086, (34) = 1.1924372702961963, (35) = 1.2129338175389266, (36) = 1.2323858241405192, (37) = 1.250955152949817, (38) = 1.268795762648431, (39) = 1.2859897085262963, (40) = 1.3026146566339305, (41) = 1.3187452976919192, (42) = 1.3344237627757964, (43) = 1.3497030677053556, (44) = 1.3646567912885574, (45) = 1.3793073520205892, (46) = 1.3936816310415334, (47) = 1.407863413800478, (48) = 1.421874455055364, (49) = 1.4357209480594326, (50) = 1.449452395680341, (51) = 1.463105461314161, (52) = 1.4766814720936998, (53) = 1.4901926102416252, (54) = 1.5036746366401115, (55) = 1.5171303716739315, (56) = 1.5305608572874234, (57) = 1.5439797161303026, (58) = 1.5573915356105725, (59) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(59, 2, {(1, 1) = .0, (1, 2) = 12.679746440910824, (2, 1) = .17295641157725689, (2, 2) = 12.567764844010108, (3, 1) = .3503687425101184, (3, 2) = 12.226735465168346, (4, 1) = .5296761539769358, (4, 2) = 11.66918438963109, (5, 1) = .7057258717088374, (5, 2) = 10.944263938053027, (6, 1) = .8727045160581467, (6, 2) = 10.12940781981862, (7, 1) = 1.0293675300859457, (7, 2) = 9.286027153268176, (8, 1) = 1.1758543805831463, (8, 2) = 8.457430747175952, (9, 1) = 1.3127197264585053, (9, 2) = 7.671647779381384, (10, 1) = 1.4404597517069202, (10, 2) = 6.945889507639623, (11, 1) = 1.5602083118771684, (11, 2) = 6.28529164237151, (12, 1) = 1.6729138904127339, (12, 2) = 5.690377357385402, (13, 1) = 1.7796245072672423, (13, 2) = 5.157666928630136, (14, 1) = 1.8817159269340347, (14, 2) = 4.680400525330798, (15, 1) = 1.9802286821314778, (15, 2) = 4.2529853303192855, (16, 1) = 2.0768238867507085, (16, 2) = 3.8675310705673276, (17, 1) = 2.173068367774721, (17, 2) = 3.5178917521371798, (18, 1) = 2.270890048473742, (18, 2) = 3.1984139900960797, (19, 1) = 2.3732345057338295, (19, 2) = 2.9030738541180425, (20, 1) = 2.4832304722226684, (20, 2) = 2.6296687269685846, (21, 1) = 2.607217350093064, (21, 2) = 2.3751465147467514, (22, 1) = 2.7547856154152726, (22, 2) = 2.144274550259719, (23, 1) = 2.939616093527199, (23, 2) = 1.961813066469559, (24, 1) = 3.162337834655888, (24, 2) = 1.8947028342219836, (25, 1) = 3.3794261487929442, (25, 2) = 1.9881175251450118, (26, 1) = 3.5551825581028798, (26, 2) = 2.1804710450280265, (27, 1) = 3.69665130286631, (27, 2) = 2.41371269610685, (28, 1) = 3.8169541259171704, (28, 2) = 2.6689707564501304, (29, 1) = 3.9243694765971413, (29, 2) = 2.9422812417392152, (30, 1) = 4.024767234931778, (30, 2) = 3.237123523436467, (31, 1) = 4.121259324516239, (31, 2) = 3.5565852574603505, (32, 1) = 4.216529428559537, (32, 2) = 3.9064929391469123, (33, 1) = 4.312457356149003, (33, 2) = 4.292723456572624, (34, 1) = 4.410412203730131, (34, 2) = 4.720850296041311, (35, 1) = 4.511931211523759, (35, 2) = 5.198255733142234, (36, 1) = 4.618090131180544, (36, 2) = 5.730648052338743, (37, 1) = 4.729865319802914, (37, 2) = 6.322579877258892, (38, 1) = 4.848375884647389, (38, 2) = 6.977624814470957, (39, 1) = 4.974387707359408, (39, 2) = 7.694131606716463, (40, 1) = 5.108607237299307, (40, 2) = 8.464739698798486, (41, 1) = 5.251606222026279, (41, 2) = 9.273734329435953, (42, 1) = 5.403405586640311, (42, 2) = 10.092710388426758, (43, 1) = 5.563684034528413, (43, 2) = 10.88121631033313, (44, 1) = 5.731806698053641, (44, 2) = 11.588598491150355, (45, 1) = 5.905935538945928, (45, 2) = 12.156134835319502, (46, 1) = 6.0836133491883215, (46, 2) = 12.530876289186718, (47, 1) = 6.262647078633722, (47, 2) = 12.678160168039451, (48, 1) = 6.439921199996006, (48, 2) = 12.587708080348033, (49, 1) = 6.612315033049238, (49, 2) = 12.2791039631489, (50, 1) = 6.777770156025198, (50, 2) = 11.79378728978802, (51, 1) = 6.934742040084409, (51, 2) = 11.183836634124887, (52, 1) = 7.0820078226929555, (52, 2) = 10.502833198232775, (53, 1) = 7.219149840143932, (53, 2) = 9.796155822525524, (54, 1) = 7.34648415514996, (54, 2) = 9.096522907410908, (55, 1) = 7.464338781940853, (55, 2) = 8.427075143053173, (56, 1) = 7.573265028033481, (56, 2) = 7.801521511288004, (57, 1) = 7.67403328840954, (57, 2) = 7.22588851589702, (58, 1) = 7.767374215067717, (58, 2) = 6.701837954309735, (59, 1) = 7.85398163397448, (59, 2) = 6.228168238234758}, datatype = float[8], order = C_order); YP := Matrix(59, 2, {(1, 1) = 12.679746440910824, (1, 2) = -.0, (2, 1) = 12.567764844010108, (2, 2) = -16.198974726732907, (3, 1) = 12.226735465168346, (3, 2) = -31.019959257656904, (4, 1) = 11.66918438963109, (4, 2) = -42.61556998298658, (5, 1) = 10.944263938053027, (5, 2) = -49.758388561480395, (6, 1) = 10.12940781981862, (6, 2) = -52.428967468700804, (7, 1) = 9.286027153268176, (7, 2) = -51.59782241363348, (8, 1) = 8.457430747175952, (8, 2) = -48.43560278055239, (9, 1) = 7.671647779381384, (9, 2) = -43.97106435599, (10, 1) = 6.945889507639623, (10, 2) = -38.985423814257885, (11, 1) = 6.28529164237151, (11, 2) = -33.97931379570107, (12, 1) = 5.690377357385402, (12, 2) = -29.25037640167964, (13, 1) = 5.157666928630136, (13, 2) = -24.944278873434246, (14, 1) = 4.680400525330798, (14, 2) = -21.099971062135964, (15, 1) = 4.2529853303192855, (15, 2) = -17.71654087434397, (16, 1) = 3.8675310705673276, (16, 2) = -14.74668416472017, (17, 1) = 3.5178917521371798, (17, 2) = -12.14221532932914, (18, 1) = 3.1984139900960797, (18, 2) = -9.851136892783515, (19, 1) = 2.9030738541180425, (19, 2) = -7.815047340827251, (20, 1) = 2.6296687269685846, (20, 2) = -5.996847641998252, (21, 1) = 2.3751465147467514, (21, 2) = -4.340048329343317, (22, 1) = 2.144274550259719, (22, 2) = -2.7978073376293024, (23, 1) = 1.961813066469559, (23, 2) = -1.320655805545758, (24, 1) = 1.8947028342219836, (24, 2) = .13039610824419176, (25, 1) = 1.9881175251450118, (25, 2) = 1.5787914768274436, (26, 1) = 2.1804710450280265, (26, 2) = 3.0488939854293733, (27, 1) = 2.41371269610685, (27, 2) = 4.590848716067672, (28, 1) = 2.6689707564501304, (28, 2) = 6.254836311594842, (29, 1) = 2.9422812417392152, (29, 2) = 8.080806575409401, (30, 1) = 3.237123523436467, (30, 2) = 10.12389716847413, (31, 1) = 3.5565852574603505, (31, 2) = 12.425685648274625, (32, 1) = 3.9064929391469123, (32, 2) = 15.042554769209124, (33, 1) = 4.292723456572624, (33, 2) = 18.027628684062055, (34, 1) = 4.720850296041311, (34, 2) = 21.42377820734578, (35, 1) = 5.198255733142234, (35, 2) = 25.27273224998513, (36, 1) = 5.730648052338743, (36, 2) = 29.574358687095494, (37, 1) = 6.322579877258892, (37, 2) = 34.27026010854073, (38, 1) = 6.977624814470957, (38, 2) = 39.21645062769387, (39, 1) = 7.694131606716463, (39, 2) = 44.11406880450838, (40, 1) = 8.464739698798486, (40, 2) = 48.471285191096335, (41, 1) = 9.273734329435953, (41, 2) = 51.56564318350719, (42, 1) = 10.092710388426758, (42, 2) = 52.45733113498812, (43, 1) = 10.88121631033313, (43, 2) = 50.12848157993341, (44, 1) = 11.588598491150355, (44, 2) = 43.73470437900077, (45, 1) = 12.156134835319502, (45, 2) = 33.00620482181849, (46, 1) = 12.530876289186718, (46, 2) = 18.57982197925053, (47, 1) = 12.678160168039451, (47, 2) = 1.9582711648426194, (48, 1) = 12.587708080348033, (48, 2) = -14.727409356953261, (49, 1) = 12.2791039631489, (49, 2) = -29.39558090446206, (50, 1) = 11.79378728978802, (50, 2) = -40.669372493045, (51, 1) = 11.183836634124887, (51, 2) = -48.03607321421275, (52, 1) = 10.502833198232775, (52, 2) = -51.72939943874306, (53, 1) = 9.796155822525524, (53, 2) = -52.452900134134595, (54, 1) = 9.096522907410908, (54, 2) = -51.04757024201683, (55, 1) = 8.427075143053173, (55, 2) = -48.28611154332724, (56, 1) = 7.801521511288004, (56, 2) = -44.78604220693657, (57, 1) = 7.22588851589702, (57, 2) = -40.986432776247455, (58, 1) = 6.701837954309735, (58, 2) = -37.176690679157566, (59, 1) = 6.228168238234758, (59, 2) = -33.53206368300497}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(59, {(1) = .0, (2) = 0.13680812436809657e-1, (3) = 0.27971051145644167e-1, (4) = 0.4296023135974479e-1, (5) = 0.5851786182852681e-1, (6) = 0.7435973609701907e-1, (7) = 0.9049963296106847e-1, (8) = .1070194360291177, (9) = .12400404122395811, (10) = .14149911102616744, (11) = .15962080170244286, (12) = .17846661163773977, (13) = .1981659899247577, (14) = .2189484735830971, (15) = .2410342106109193, (16) = .2648591917099649, (17) = .29096256417430655, (18) = .3201405168424954, (19) = .3537506263213988, (20) = .39359925436482457, (21) = .44328213017143814, (22) = .508829664760916, (23) = .5993459008179163, (24) = .7156934169226271, (25) = .8282966673508598, (26) = .9130409367878883, (27) = .9748398788122531, (28) = 1.0223008015052495, (29) = 1.0606661846616776, (30) = 1.0932188074676115, (31) = 1.1216707736458404, (32) = 1.1472399927159063, (33) = 1.1706725764522086, (34) = 1.1924372702961963, (35) = 1.2129338175389266, (36) = 1.2323858241405192, (37) = 1.250955152949817, (38) = 1.268795762648431, (39) = 1.2859897085262963, (40) = 1.3026146566339305, (41) = 1.3187452976919192, (42) = 1.3344237627757964, (43) = 1.3497030677053556, (44) = 1.3646567912885574, (45) = 1.3793073520205892, (46) = 1.3936816310415334, (47) = 1.407863413800478, (48) = 1.421874455055364, (49) = 1.4357209480594326, (50) = 1.449452395680341, (51) = 1.463105461314161, (52) = 1.4766814720936998, (53) = 1.4901926102416252, (54) = 1.5036746366401115, (55) = 1.5171303716739315, (56) = 1.5305608572874234, (57) = 1.5439797161303026, (58) = 1.5573915356105725, (59) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(59, 2, {(1, 1) = .0, (1, 2) = -0.2040971830865995e-12, (2, 1) = -0.2515085127931946e-15, (2, 2) = -0.2182240868305703e-12, (3, 1) = -0.3034400215680231e-14, (3, 2) = -0.24629572182126297e-12, (4, 1) = -0.70293294530316345e-14, (4, 2) = -0.27355561891671303e-12, (5, 1) = -0.11642509457266535e-13, (5, 2) = -0.2921725988376451e-12, (6, 1) = -0.14891698174826332e-13, (6, 2) = -0.27332966679506926e-12, (7, 1) = -0.2244158402930937e-13, (7, 2) = -0.21555448086863683e-12, (8, 1) = -0.23830070405255755e-13, (8, 2) = -0.14452163468605652e-12, (9, 1) = -0.28784241689197694e-13, (9, 2) = -0.8936429016769738e-13, (10, 1) = -0.2870763992607896e-13, (10, 2) = -0.7015153701410498e-13, (11, 1) = -0.29805232120338505e-13, (11, 2) = -0.20455226198960395e-13, (12, 1) = -0.3071869736898188e-13, (12, 2) = -0.949608712395301e-14, (13, 1) = -0.33085750279284234e-13, (13, 2) = -0.9853725522439467e-14, (14, 1) = -0.3024025650769599e-13, (14, 2) = 0.55343272382320285e-14, (15, 1) = -0.2989762453456761e-13, (15, 2) = 0.9484091577415516e-14, (16, 1) = -0.3126420836975352e-13, (16, 2) = 0.10993914919372806e-15, (17, 1) = -0.2904952409517121e-13, (17, 2) = 0.7741875869166384e-14, (18, 1) = -0.26389405390354672e-13, (18, 2) = -0.39736138956257974e-14, (19, 1) = -0.3066642976676729e-13, (19, 2) = -0.8841791351636723e-14, (20, 1) = -0.3381874764737987e-13, (20, 2) = -0.10474029579197825e-13, (21, 1) = -0.3076950762715862e-13, (21, 2) = -0.14630357321113296e-13, (22, 1) = -0.25742990878046857e-13, (22, 2) = -0.32013707096004785e-13, (23, 1) = -0.26034531106847924e-14, (23, 2) = -0.17882398235769704e-13, (24, 1) = -0.24294629294277565e-13, (24, 2) = 0.6033751098446911e-13, (25, 1) = -0.10212795741481291e-13, (25, 2) = -0.5409909790245115e-14, (26, 1) = 0.65234775286404595e-14, (26, 2) = 0.10506976038853487e-14, (27, 1) = 0.15145921627767115e-13, (27, 2) = 0.11062582034309036e-13, (28, 1) = 0.17727034515154246e-13, (28, 2) = 0.1392374003706416e-13, (29, 1) = 0.20190304444405087e-13, (29, 2) = 0.24077662756883795e-13, (30, 1) = 0.15567756982704383e-13, (30, 2) = 0.3220356728033926e-13, (31, 1) = 0.2049236289955625e-13, (31, 2) = 0.4582315826902611e-13, (32, 1) = 0.2093779309445625e-13, (32, 2) = 0.50847503530540665e-13, (33, 1) = 0.34656722670579105e-13, (33, 2) = 0.65791408284367e-13, (34, 1) = 0.1414627403891834e-13, (34, 2) = 0.7434970442839986e-13, (35, 1) = 0.1430197385087178e-13, (35, 2) = 0.807995248637132e-13, (36, 1) = 0.2674246752632126e-13, (36, 2) = 0.8238316242997303e-13, (37, 1) = 0.3749521089102351e-13, (37, 2) = 0.9551972331079688e-13, (38, 1) = 0.43798173157104124e-13, (38, 2) = 0.8882243056529444e-13, (39, 1) = 0.34765247192333954e-13, (39, 2) = 0.7174268116999525e-13, (40, 1) = 0.3492718539460031e-13, (40, 2) = 0.25074509668472094e-13, (41, 1) = 0.3268343950774218e-13, (41, 2) = -0.6182099171300564e-13, (42, 1) = 0.7615127467050659e-14, (42, 2) = -0.1854246486047524e-12, (43, 1) = 0.20280258788197414e-13, (43, 2) = -0.2612655964015281e-12, (44, 1) = 0.5268977580072408e-14, (44, 2) = -0.3058713868413555e-12, (45, 1) = 0.8243135377351008e-14, (45, 2) = -0.3000324869932838e-12, (46, 1) = 0.9651203882669847e-14, (46, 2) = -0.23958986839330786e-12, (47, 1) = -0.15172603499122148e-13, (47, 2) = -0.14114952424844551e-12, (48, 1) = 0.6415293117671377e-14, (48, 2) = -0.3574626525598719e-13, (49, 1) = 0.11771473310921094e-13, (49, 2) = 0.12758808412240534e-13, (50, 1) = 0.350361321759921e-14, (50, 2) = 0.18960251337597162e-13, (51, 1) = -0.402567915071639e-14, (51, 2) = -0.42598728779319296e-13, (52, 1) = 0.14493706543893737e-13, (52, 2) = -0.4380811741884113e-13, (53, 1) = 0.19531436506149837e-13, (53, 2) = -0.8490114617635297e-13, (54, 1) = 0.4784180312117899e-14, (54, 2) = -0.10238102255674353e-12, (55, 1) = 0.13186433513432542e-13, (55, 2) = -0.12142340951879697e-12, (56, 1) = 0.68589188913736246e-14, (56, 2) = -0.10734720880680753e-12, (57, 1) = 0.4428481028761533e-14, (57, 2) = -0.9866461770792047e-13, (58, 1) = 0.7387414481358141e-14, (58, 2) = -0.9166293971010342e-13, (59, 1) = .0, (59, 2) = -0.731647525953621e-13}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[59] elif outpoint = "order" then return 10 elif outpoint = "error" then return HFloat(3.058713868413555e-13) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 59, [v(u), diff(v(u), u)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[59] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[59] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(59, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(59, 2, X, Y, outpoint, yout, L, V) end if; [u = outpoint, seq('[v(u), diff(v(u), u)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[59] elif outpoint = "order" then return 10 elif outpoint = "error" then return HFloat(3.058713868413555e-13) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 59, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[59] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[59] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(59, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(59, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 36893628625353841228, (2) = 36893628625353841404, (3) = 36893628625353841580}), (3) = [u, v(u), diff(v(u), u)], (4) = 0}); _solnproc := _dat[1]; if member(u, ["last", 'last']) then _res := _solnproc("last"); if type(_res, 'list') then return _res[1] end if elif type(u, `=`) and member(lhs(u), ["initial", 'initial']) then if type(rhs(u), 'list') then _res := _solnproc("initial" = [0, op(rhs(u))]) else _res := _solnproc("initial" = [1, rhs(u)]) end if; if type(_res, 'list') then return _res[1] end if elif u = "sysvars" then return _dat[3] end if; u end proc, v = proc (u) local res, data, solnproc, v, outpoint; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](u) else outpoint := evalf(u) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(59, {(1) = .0, (2) = 0.13680812436809657e-1, (3) = 0.27971051145644167e-1, (4) = 0.4296023135974479e-1, (5) = 0.5851786182852681e-1, (6) = 0.7435973609701907e-1, (7) = 0.9049963296106847e-1, (8) = .1070194360291177, (9) = .12400404122395811, (10) = .14149911102616744, (11) = .15962080170244286, (12) = .17846661163773977, (13) = .1981659899247577, (14) = .2189484735830971, (15) = .2410342106109193, (16) = .2648591917099649, (17) = .29096256417430655, (18) = .3201405168424954, (19) = .3537506263213988, (20) = .39359925436482457, (21) = .44328213017143814, (22) = .508829664760916, (23) = .5993459008179163, (24) = .7156934169226271, (25) = .8282966673508598, (26) = .9130409367878883, (27) = .9748398788122531, (28) = 1.0223008015052495, (29) = 1.0606661846616776, (30) = 1.0932188074676115, (31) = 1.1216707736458404, (32) = 1.1472399927159063, (33) = 1.1706725764522086, (34) = 1.1924372702961963, (35) = 1.2129338175389266, (36) = 1.2323858241405192, (37) = 1.250955152949817, (38) = 1.268795762648431, (39) = 1.2859897085262963, (40) = 1.3026146566339305, (41) = 1.3187452976919192, (42) = 1.3344237627757964, (43) = 1.3497030677053556, (44) = 1.3646567912885574, (45) = 1.3793073520205892, (46) = 1.3936816310415334, (47) = 1.407863413800478, (48) = 1.421874455055364, (49) = 1.4357209480594326, (50) = 1.449452395680341, (51) = 1.463105461314161, (52) = 1.4766814720936998, (53) = 1.4901926102416252, (54) = 1.5036746366401115, (55) = 1.5171303716739315, (56) = 1.5305608572874234, (57) = 1.5439797161303026, (58) = 1.5573915356105725, (59) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(59, 2, {(1, 1) = .0, (1, 2) = 12.679746440910824, (2, 1) = .17295641157725689, (2, 2) = 12.567764844010108, (3, 1) = .3503687425101184, (3, 2) = 12.226735465168346, (4, 1) = .5296761539769358, (4, 2) = 11.66918438963109, (5, 1) = .7057258717088374, (5, 2) = 10.944263938053027, (6, 1) = .8727045160581467, (6, 2) = 10.12940781981862, (7, 1) = 1.0293675300859457, (7, 2) = 9.286027153268176, (8, 1) = 1.1758543805831463, (8, 2) = 8.457430747175952, (9, 1) = 1.3127197264585053, (9, 2) = 7.671647779381384, (10, 1) = 1.4404597517069202, (10, 2) = 6.945889507639623, (11, 1) = 1.5602083118771684, (11, 2) = 6.28529164237151, (12, 1) = 1.6729138904127339, (12, 2) = 5.690377357385402, (13, 1) = 1.7796245072672423, (13, 2) = 5.157666928630136, (14, 1) = 1.8817159269340347, (14, 2) = 4.680400525330798, (15, 1) = 1.9802286821314778, (15, 2) = 4.2529853303192855, (16, 1) = 2.0768238867507085, (16, 2) = 3.8675310705673276, (17, 1) = 2.173068367774721, (17, 2) = 3.5178917521371798, (18, 1) = 2.270890048473742, (18, 2) = 3.1984139900960797, (19, 1) = 2.3732345057338295, (19, 2) = 2.9030738541180425, (20, 1) = 2.4832304722226684, (20, 2) = 2.6296687269685846, (21, 1) = 2.607217350093064, (21, 2) = 2.3751465147467514, (22, 1) = 2.7547856154152726, (22, 2) = 2.144274550259719, (23, 1) = 2.939616093527199, (23, 2) = 1.961813066469559, (24, 1) = 3.162337834655888, (24, 2) = 1.8947028342219836, (25, 1) = 3.3794261487929442, (25, 2) = 1.9881175251450118, (26, 1) = 3.5551825581028798, (26, 2) = 2.1804710450280265, (27, 1) = 3.69665130286631, (27, 2) = 2.41371269610685, (28, 1) = 3.8169541259171704, (28, 2) = 2.6689707564501304, (29, 1) = 3.9243694765971413, (29, 2) = 2.9422812417392152, (30, 1) = 4.024767234931778, (30, 2) = 3.237123523436467, (31, 1) = 4.121259324516239, (31, 2) = 3.5565852574603505, (32, 1) = 4.216529428559537, (32, 2) = 3.9064929391469123, (33, 1) = 4.312457356149003, (33, 2) = 4.292723456572624, (34, 1) = 4.410412203730131, (34, 2) = 4.720850296041311, (35, 1) = 4.511931211523759, (35, 2) = 5.198255733142234, (36, 1) = 4.618090131180544, (36, 2) = 5.730648052338743, (37, 1) = 4.729865319802914, (37, 2) = 6.322579877258892, (38, 1) = 4.848375884647389, (38, 2) = 6.977624814470957, (39, 1) = 4.974387707359408, (39, 2) = 7.694131606716463, (40, 1) = 5.108607237299307, (40, 2) = 8.464739698798486, (41, 1) = 5.251606222026279, (41, 2) = 9.273734329435953, (42, 1) = 5.403405586640311, (42, 2) = 10.092710388426758, (43, 1) = 5.563684034528413, (43, 2) = 10.88121631033313, (44, 1) = 5.731806698053641, (44, 2) = 11.588598491150355, (45, 1) = 5.905935538945928, (45, 2) = 12.156134835319502, (46, 1) = 6.0836133491883215, (46, 2) = 12.530876289186718, (47, 1) = 6.262647078633722, (47, 2) = 12.678160168039451, (48, 1) = 6.439921199996006, (48, 2) = 12.587708080348033, (49, 1) = 6.612315033049238, (49, 2) = 12.2791039631489, (50, 1) = 6.777770156025198, (50, 2) = 11.79378728978802, (51, 1) = 6.934742040084409, (51, 2) = 11.183836634124887, (52, 1) = 7.0820078226929555, (52, 2) = 10.502833198232775, (53, 1) = 7.219149840143932, (53, 2) = 9.796155822525524, (54, 1) = 7.34648415514996, (54, 2) = 9.096522907410908, (55, 1) = 7.464338781940853, (55, 2) = 8.427075143053173, (56, 1) = 7.573265028033481, (56, 2) = 7.801521511288004, (57, 1) = 7.67403328840954, (57, 2) = 7.22588851589702, (58, 1) = 7.767374215067717, (58, 2) = 6.701837954309735, (59, 1) = 7.85398163397448, (59, 2) = 6.228168238234758}, datatype = float[8], order = C_order); YP := Matrix(59, 2, {(1, 1) = 12.679746440910824, (1, 2) = -.0, (2, 1) = 12.567764844010108, (2, 2) = -16.198974726732907, (3, 1) = 12.226735465168346, (3, 2) = -31.019959257656904, (4, 1) = 11.66918438963109, (4, 2) = -42.61556998298658, (5, 1) = 10.944263938053027, (5, 2) = -49.758388561480395, (6, 1) = 10.12940781981862, (6, 2) = -52.428967468700804, (7, 1) = 9.286027153268176, (7, 2) = -51.59782241363348, (8, 1) = 8.457430747175952, (8, 2) = -48.43560278055239, (9, 1) = 7.671647779381384, (9, 2) = -43.97106435599, (10, 1) = 6.945889507639623, (10, 2) = -38.985423814257885, (11, 1) = 6.28529164237151, (11, 2) = -33.97931379570107, (12, 1) = 5.690377357385402, (12, 2) = -29.25037640167964, (13, 1) = 5.157666928630136, (13, 2) = -24.944278873434246, (14, 1) = 4.680400525330798, (14, 2) = -21.099971062135964, (15, 1) = 4.2529853303192855, (15, 2) = -17.71654087434397, (16, 1) = 3.8675310705673276, (16, 2) = -14.74668416472017, (17, 1) = 3.5178917521371798, (17, 2) = -12.14221532932914, (18, 1) = 3.1984139900960797, (18, 2) = -9.851136892783515, (19, 1) = 2.9030738541180425, (19, 2) = -7.815047340827251, (20, 1) = 2.6296687269685846, (20, 2) = -5.996847641998252, (21, 1) = 2.3751465147467514, (21, 2) = -4.340048329343317, (22, 1) = 2.144274550259719, (22, 2) = -2.7978073376293024, (23, 1) = 1.961813066469559, (23, 2) = -1.320655805545758, (24, 1) = 1.8947028342219836, (24, 2) = .13039610824419176, (25, 1) = 1.9881175251450118, (25, 2) = 1.5787914768274436, (26, 1) = 2.1804710450280265, (26, 2) = 3.0488939854293733, (27, 1) = 2.41371269610685, (27, 2) = 4.590848716067672, (28, 1) = 2.6689707564501304, (28, 2) = 6.254836311594842, (29, 1) = 2.9422812417392152, (29, 2) = 8.080806575409401, (30, 1) = 3.237123523436467, (30, 2) = 10.12389716847413, (31, 1) = 3.5565852574603505, (31, 2) = 12.425685648274625, (32, 1) = 3.9064929391469123, (32, 2) = 15.042554769209124, (33, 1) = 4.292723456572624, (33, 2) = 18.027628684062055, (34, 1) = 4.720850296041311, (34, 2) = 21.42377820734578, (35, 1) = 5.198255733142234, (35, 2) = 25.27273224998513, (36, 1) = 5.730648052338743, (36, 2) = 29.574358687095494, (37, 1) = 6.322579877258892, (37, 2) = 34.27026010854073, (38, 1) = 6.977624814470957, (38, 2) = 39.21645062769387, (39, 1) = 7.694131606716463, (39, 2) = 44.11406880450838, (40, 1) = 8.464739698798486, (40, 2) = 48.471285191096335, (41, 1) = 9.273734329435953, (41, 2) = 51.56564318350719, (42, 1) = 10.092710388426758, (42, 2) = 52.45733113498812, (43, 1) = 10.88121631033313, (43, 2) = 50.12848157993341, (44, 1) = 11.588598491150355, (44, 2) = 43.73470437900077, (45, 1) = 12.156134835319502, (45, 2) = 33.00620482181849, (46, 1) = 12.530876289186718, (46, 2) = 18.57982197925053, (47, 1) = 12.678160168039451, (47, 2) = 1.9582711648426194, (48, 1) = 12.587708080348033, (48, 2) = -14.727409356953261, (49, 1) = 12.2791039631489, (49, 2) = -29.39558090446206, (50, 1) = 11.79378728978802, (50, 2) = -40.669372493045, (51, 1) = 11.183836634124887, (51, 2) = -48.03607321421275, (52, 1) = 10.502833198232775, (52, 2) = -51.72939943874306, (53, 1) = 9.796155822525524, (53, 2) = -52.452900134134595, (54, 1) = 9.096522907410908, (54, 2) = -51.04757024201683, (55, 1) = 8.427075143053173, (55, 2) = -48.28611154332724, (56, 1) = 7.801521511288004, (56, 2) = -44.78604220693657, (57, 1) = 7.22588851589702, (57, 2) = -40.986432776247455, (58, 1) = 6.701837954309735, (58, 2) = -37.176690679157566, (59, 1) = 6.228168238234758, (59, 2) = -33.53206368300497}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(59, {(1) = .0, (2) = 0.13680812436809657e-1, (3) = 0.27971051145644167e-1, (4) = 0.4296023135974479e-1, (5) = 0.5851786182852681e-1, (6) = 0.7435973609701907e-1, (7) = 0.9049963296106847e-1, (8) = .1070194360291177, (9) = .12400404122395811, (10) = .14149911102616744, (11) = .15962080170244286, (12) = .17846661163773977, (13) = .1981659899247577, (14) = .2189484735830971, (15) = .2410342106109193, (16) = .2648591917099649, (17) = .29096256417430655, (18) = .3201405168424954, (19) = .3537506263213988, (20) = .39359925436482457, (21) = .44328213017143814, (22) = .508829664760916, (23) = .5993459008179163, (24) = .7156934169226271, (25) = .8282966673508598, (26) = .9130409367878883, (27) = .9748398788122531, (28) = 1.0223008015052495, (29) = 1.0606661846616776, (30) = 1.0932188074676115, (31) = 1.1216707736458404, (32) = 1.1472399927159063, (33) = 1.1706725764522086, (34) = 1.1924372702961963, (35) = 1.2129338175389266, (36) = 1.2323858241405192, (37) = 1.250955152949817, (38) = 1.268795762648431, (39) = 1.2859897085262963, (40) = 1.3026146566339305, (41) = 1.3187452976919192, (42) = 1.3344237627757964, (43) = 1.3497030677053556, (44) = 1.3646567912885574, (45) = 1.3793073520205892, (46) = 1.3936816310415334, (47) = 1.407863413800478, (48) = 1.421874455055364, (49) = 1.4357209480594326, (50) = 1.449452395680341, (51) = 1.463105461314161, (52) = 1.4766814720936998, (53) = 1.4901926102416252, (54) = 1.5036746366401115, (55) = 1.5171303716739315, (56) = 1.5305608572874234, (57) = 1.5439797161303026, (58) = 1.5573915356105725, (59) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(59, 2, {(1, 1) = .0, (1, 2) = -0.2040971830865995e-12, (2, 1) = -0.2515085127931946e-15, (2, 2) = -0.2182240868305703e-12, (3, 1) = -0.3034400215680231e-14, (3, 2) = -0.24629572182126297e-12, (4, 1) = -0.70293294530316345e-14, (4, 2) = -0.27355561891671303e-12, (5, 1) = -0.11642509457266535e-13, (5, 2) = -0.2921725988376451e-12, (6, 1) = -0.14891698174826332e-13, (6, 2) = -0.27332966679506926e-12, (7, 1) = -0.2244158402930937e-13, (7, 2) = -0.21555448086863683e-12, (8, 1) = -0.23830070405255755e-13, (8, 2) = -0.14452163468605652e-12, (9, 1) = -0.28784241689197694e-13, (9, 2) = -0.8936429016769738e-13, (10, 1) = -0.2870763992607896e-13, (10, 2) = -0.7015153701410498e-13, (11, 1) = -0.29805232120338505e-13, (11, 2) = -0.20455226198960395e-13, (12, 1) = -0.3071869736898188e-13, (12, 2) = -0.949608712395301e-14, (13, 1) = -0.33085750279284234e-13, (13, 2) = -0.9853725522439467e-14, (14, 1) = -0.3024025650769599e-13, (14, 2) = 0.55343272382320285e-14, (15, 1) = -0.2989762453456761e-13, (15, 2) = 0.9484091577415516e-14, (16, 1) = -0.3126420836975352e-13, (16, 2) = 0.10993914919372806e-15, (17, 1) = -0.2904952409517121e-13, (17, 2) = 0.7741875869166384e-14, (18, 1) = -0.26389405390354672e-13, (18, 2) = -0.39736138956257974e-14, (19, 1) = -0.3066642976676729e-13, (19, 2) = -0.8841791351636723e-14, (20, 1) = -0.3381874764737987e-13, (20, 2) = -0.10474029579197825e-13, (21, 1) = -0.3076950762715862e-13, (21, 2) = -0.14630357321113296e-13, (22, 1) = -0.25742990878046857e-13, (22, 2) = -0.32013707096004785e-13, (23, 1) = -0.26034531106847924e-14, (23, 2) = -0.17882398235769704e-13, (24, 1) = -0.24294629294277565e-13, (24, 2) = 0.6033751098446911e-13, (25, 1) = -0.10212795741481291e-13, (25, 2) = -0.5409909790245115e-14, (26, 1) = 0.65234775286404595e-14, (26, 2) = 0.10506976038853487e-14, (27, 1) = 0.15145921627767115e-13, (27, 2) = 0.11062582034309036e-13, (28, 1) = 0.17727034515154246e-13, (28, 2) = 0.1392374003706416e-13, (29, 1) = 0.20190304444405087e-13, (29, 2) = 0.24077662756883795e-13, (30, 1) = 0.15567756982704383e-13, (30, 2) = 0.3220356728033926e-13, (31, 1) = 0.2049236289955625e-13, (31, 2) = 0.4582315826902611e-13, (32, 1) = 0.2093779309445625e-13, (32, 2) = 0.50847503530540665e-13, (33, 1) = 0.34656722670579105e-13, (33, 2) = 0.65791408284367e-13, (34, 1) = 0.1414627403891834e-13, (34, 2) = 0.7434970442839986e-13, (35, 1) = 0.1430197385087178e-13, (35, 2) = 0.807995248637132e-13, (36, 1) = 0.2674246752632126e-13, (36, 2) = 0.8238316242997303e-13, (37, 1) = 0.3749521089102351e-13, (37, 2) = 0.9551972331079688e-13, (38, 1) = 0.43798173157104124e-13, (38, 2) = 0.8882243056529444e-13, (39, 1) = 0.34765247192333954e-13, (39, 2) = 0.7174268116999525e-13, (40, 1) = 0.3492718539460031e-13, (40, 2) = 0.25074509668472094e-13, (41, 1) = 0.3268343950774218e-13, (41, 2) = -0.6182099171300564e-13, (42, 1) = 0.7615127467050659e-14, (42, 2) = -0.1854246486047524e-12, (43, 1) = 0.20280258788197414e-13, (43, 2) = -0.2612655964015281e-12, (44, 1) = 0.5268977580072408e-14, (44, 2) = -0.3058713868413555e-12, (45, 1) = 0.8243135377351008e-14, (45, 2) = -0.3000324869932838e-12, (46, 1) = 0.9651203882669847e-14, (46, 2) = -0.23958986839330786e-12, (47, 1) = -0.15172603499122148e-13, (47, 2) = -0.14114952424844551e-12, (48, 1) = 0.6415293117671377e-14, (48, 2) = -0.3574626525598719e-13, (49, 1) = 0.11771473310921094e-13, (49, 2) = 0.12758808412240534e-13, (50, 1) = 0.350361321759921e-14, (50, 2) = 0.18960251337597162e-13, (51, 1) = -0.402567915071639e-14, (51, 2) = -0.42598728779319296e-13, (52, 1) = 0.14493706543893737e-13, (52, 2) = -0.4380811741884113e-13, (53, 1) = 0.19531436506149837e-13, (53, 2) = -0.8490114617635297e-13, (54, 1) = 0.4784180312117899e-14, (54, 2) = -0.10238102255674353e-12, (55, 1) = 0.13186433513432542e-13, (55, 2) = -0.12142340951879697e-12, (56, 1) = 0.68589188913736246e-14, (56, 2) = -0.10734720880680753e-12, (57, 1) = 0.4428481028761533e-14, (57, 2) = -0.9866461770792047e-13, (58, 1) = 0.7387414481358141e-14, (58, 2) = -0.9166293971010342e-13, (59, 1) = .0, (59, 2) = -0.731647525953621e-13}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[59] elif outpoint = "order" then return 10 elif outpoint = "error" then return HFloat(3.058713868413555e-13) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 59, [v(u), diff(v(u), u)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[59] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[59] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(59, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(59, 2, X, Y, outpoint, yout, L, V) end if; [u = outpoint, seq('[v(u), diff(v(u), u)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[59] elif outpoint = "order" then return 10 elif outpoint = "error" then return HFloat(3.058713868413555e-13) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 59, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[59] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[59] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(59, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(59, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 36893628625353841228, (2) = 36893628625353841404, (3) = 36893628625353841580}), (3) = [u, v(u), diff(v(u), u)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(u) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(u) else v := pointto(data[2][2]); return ('v')(u) end if end if; try res := solnproc(outpoint); res[2] catch: error  end try end proc, D(v) = proc (u) local res, data, solnproc, `D(v)`, outpoint; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](u) else outpoint := evalf(u) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(59, {(1) = .0, (2) = 0.13680812436809657e-1, (3) = 0.27971051145644167e-1, (4) = 0.4296023135974479e-1, (5) = 0.5851786182852681e-1, (6) = 0.7435973609701907e-1, (7) = 0.9049963296106847e-1, (8) = .1070194360291177, (9) = .12400404122395811, (10) = .14149911102616744, (11) = .15962080170244286, (12) = .17846661163773977, (13) = .1981659899247577, (14) = .2189484735830971, (15) = .2410342106109193, (16) = .2648591917099649, (17) = .29096256417430655, (18) = .3201405168424954, (19) = .3537506263213988, (20) = .39359925436482457, (21) = .44328213017143814, (22) = .508829664760916, (23) = .5993459008179163, (24) = .7156934169226271, (25) = .8282966673508598, (26) = .9130409367878883, (27) = .9748398788122531, (28) = 1.0223008015052495, (29) = 1.0606661846616776, (30) = 1.0932188074676115, (31) = 1.1216707736458404, (32) = 1.1472399927159063, (33) = 1.1706725764522086, (34) = 1.1924372702961963, (35) = 1.2129338175389266, (36) = 1.2323858241405192, (37) = 1.250955152949817, (38) = 1.268795762648431, (39) = 1.2859897085262963, (40) = 1.3026146566339305, (41) = 1.3187452976919192, (42) = 1.3344237627757964, (43) = 1.3497030677053556, (44) = 1.3646567912885574, (45) = 1.3793073520205892, (46) = 1.3936816310415334, (47) = 1.407863413800478, (48) = 1.421874455055364, (49) = 1.4357209480594326, (50) = 1.449452395680341, (51) = 1.463105461314161, (52) = 1.4766814720936998, (53) = 1.4901926102416252, (54) = 1.5036746366401115, (55) = 1.5171303716739315, (56) = 1.5305608572874234, (57) = 1.5439797161303026, (58) = 1.5573915356105725, (59) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(59, 2, {(1, 1) = .0, (1, 2) = 12.679746440910824, (2, 1) = .17295641157725689, (2, 2) = 12.567764844010108, (3, 1) = .3503687425101184, (3, 2) = 12.226735465168346, (4, 1) = .5296761539769358, (4, 2) = 11.66918438963109, (5, 1) = .7057258717088374, (5, 2) = 10.944263938053027, (6, 1) = .8727045160581467, (6, 2) = 10.12940781981862, (7, 1) = 1.0293675300859457, (7, 2) = 9.286027153268176, (8, 1) = 1.1758543805831463, (8, 2) = 8.457430747175952, (9, 1) = 1.3127197264585053, (9, 2) = 7.671647779381384, (10, 1) = 1.4404597517069202, (10, 2) = 6.945889507639623, (11, 1) = 1.5602083118771684, (11, 2) = 6.28529164237151, (12, 1) = 1.6729138904127339, (12, 2) = 5.690377357385402, (13, 1) = 1.7796245072672423, (13, 2) = 5.157666928630136, (14, 1) = 1.8817159269340347, (14, 2) = 4.680400525330798, (15, 1) = 1.9802286821314778, (15, 2) = 4.2529853303192855, (16, 1) = 2.0768238867507085, (16, 2) = 3.8675310705673276, (17, 1) = 2.173068367774721, (17, 2) = 3.5178917521371798, (18, 1) = 2.270890048473742, (18, 2) = 3.1984139900960797, (19, 1) = 2.3732345057338295, (19, 2) = 2.9030738541180425, (20, 1) = 2.4832304722226684, (20, 2) = 2.6296687269685846, (21, 1) = 2.607217350093064, (21, 2) = 2.3751465147467514, (22, 1) = 2.7547856154152726, (22, 2) = 2.144274550259719, (23, 1) = 2.939616093527199, (23, 2) = 1.961813066469559, (24, 1) = 3.162337834655888, (24, 2) = 1.8947028342219836, (25, 1) = 3.3794261487929442, (25, 2) = 1.9881175251450118, (26, 1) = 3.5551825581028798, (26, 2) = 2.1804710450280265, (27, 1) = 3.69665130286631, (27, 2) = 2.41371269610685, (28, 1) = 3.8169541259171704, (28, 2) = 2.6689707564501304, (29, 1) = 3.9243694765971413, (29, 2) = 2.9422812417392152, (30, 1) = 4.024767234931778, (30, 2) = 3.237123523436467, (31, 1) = 4.121259324516239, (31, 2) = 3.5565852574603505, (32, 1) = 4.216529428559537, (32, 2) = 3.9064929391469123, (33, 1) = 4.312457356149003, (33, 2) = 4.292723456572624, (34, 1) = 4.410412203730131, (34, 2) = 4.720850296041311, (35, 1) = 4.511931211523759, (35, 2) = 5.198255733142234, (36, 1) = 4.618090131180544, (36, 2) = 5.730648052338743, (37, 1) = 4.729865319802914, (37, 2) = 6.322579877258892, (38, 1) = 4.848375884647389, (38, 2) = 6.977624814470957, (39, 1) = 4.974387707359408, (39, 2) = 7.694131606716463, (40, 1) = 5.108607237299307, (40, 2) = 8.464739698798486, (41, 1) = 5.251606222026279, (41, 2) = 9.273734329435953, (42, 1) = 5.403405586640311, (42, 2) = 10.092710388426758, (43, 1) = 5.563684034528413, (43, 2) = 10.88121631033313, (44, 1) = 5.731806698053641, (44, 2) = 11.588598491150355, (45, 1) = 5.905935538945928, (45, 2) = 12.156134835319502, (46, 1) = 6.0836133491883215, (46, 2) = 12.530876289186718, (47, 1) = 6.262647078633722, (47, 2) = 12.678160168039451, (48, 1) = 6.439921199996006, (48, 2) = 12.587708080348033, (49, 1) = 6.612315033049238, (49, 2) = 12.2791039631489, (50, 1) = 6.777770156025198, (50, 2) = 11.79378728978802, (51, 1) = 6.934742040084409, (51, 2) = 11.183836634124887, (52, 1) = 7.0820078226929555, (52, 2) = 10.502833198232775, (53, 1) = 7.219149840143932, (53, 2) = 9.796155822525524, (54, 1) = 7.34648415514996, (54, 2) = 9.096522907410908, (55, 1) = 7.464338781940853, (55, 2) = 8.427075143053173, (56, 1) = 7.573265028033481, (56, 2) = 7.801521511288004, (57, 1) = 7.67403328840954, (57, 2) = 7.22588851589702, (58, 1) = 7.767374215067717, (58, 2) = 6.701837954309735, (59, 1) = 7.85398163397448, (59, 2) = 6.228168238234758}, datatype = float[8], order = C_order); YP := Matrix(59, 2, {(1, 1) = 12.679746440910824, (1, 2) = -.0, (2, 1) = 12.567764844010108, (2, 2) = -16.198974726732907, (3, 1) = 12.226735465168346, (3, 2) = -31.019959257656904, (4, 1) = 11.66918438963109, (4, 2) = -42.61556998298658, (5, 1) = 10.944263938053027, (5, 2) = -49.758388561480395, (6, 1) = 10.12940781981862, (6, 2) = -52.428967468700804, (7, 1) = 9.286027153268176, (7, 2) = -51.59782241363348, (8, 1) = 8.457430747175952, (8, 2) = -48.43560278055239, (9, 1) = 7.671647779381384, (9, 2) = -43.97106435599, (10, 1) = 6.945889507639623, (10, 2) = -38.985423814257885, (11, 1) = 6.28529164237151, (11, 2) = -33.97931379570107, (12, 1) = 5.690377357385402, (12, 2) = -29.25037640167964, (13, 1) = 5.157666928630136, (13, 2) = -24.944278873434246, (14, 1) = 4.680400525330798, (14, 2) = -21.099971062135964, (15, 1) = 4.2529853303192855, (15, 2) = -17.71654087434397, (16, 1) = 3.8675310705673276, (16, 2) = -14.74668416472017, (17, 1) = 3.5178917521371798, (17, 2) = -12.14221532932914, (18, 1) = 3.1984139900960797, (18, 2) = -9.851136892783515, (19, 1) = 2.9030738541180425, (19, 2) = -7.815047340827251, (20, 1) = 2.6296687269685846, (20, 2) = -5.996847641998252, (21, 1) = 2.3751465147467514, (21, 2) = -4.340048329343317, (22, 1) = 2.144274550259719, (22, 2) = -2.7978073376293024, (23, 1) = 1.961813066469559, (23, 2) = -1.320655805545758, (24, 1) = 1.8947028342219836, (24, 2) = .13039610824419176, (25, 1) = 1.9881175251450118, (25, 2) = 1.5787914768274436, (26, 1) = 2.1804710450280265, (26, 2) = 3.0488939854293733, (27, 1) = 2.41371269610685, (27, 2) = 4.590848716067672, (28, 1) = 2.6689707564501304, (28, 2) = 6.254836311594842, (29, 1) = 2.9422812417392152, (29, 2) = 8.080806575409401, (30, 1) = 3.237123523436467, (30, 2) = 10.12389716847413, (31, 1) = 3.5565852574603505, (31, 2) = 12.425685648274625, (32, 1) = 3.9064929391469123, (32, 2) = 15.042554769209124, (33, 1) = 4.292723456572624, (33, 2) = 18.027628684062055, (34, 1) = 4.720850296041311, (34, 2) = 21.42377820734578, (35, 1) = 5.198255733142234, (35, 2) = 25.27273224998513, (36, 1) = 5.730648052338743, (36, 2) = 29.574358687095494, (37, 1) = 6.322579877258892, (37, 2) = 34.27026010854073, (38, 1) = 6.977624814470957, (38, 2) = 39.21645062769387, (39, 1) = 7.694131606716463, (39, 2) = 44.11406880450838, (40, 1) = 8.464739698798486, (40, 2) = 48.471285191096335, (41, 1) = 9.273734329435953, (41, 2) = 51.56564318350719, (42, 1) = 10.092710388426758, (42, 2) = 52.45733113498812, (43, 1) = 10.88121631033313, (43, 2) = 50.12848157993341, (44, 1) = 11.588598491150355, (44, 2) = 43.73470437900077, (45, 1) = 12.156134835319502, (45, 2) = 33.00620482181849, (46, 1) = 12.530876289186718, (46, 2) = 18.57982197925053, (47, 1) = 12.678160168039451, (47, 2) = 1.9582711648426194, (48, 1) = 12.587708080348033, (48, 2) = -14.727409356953261, (49, 1) = 12.2791039631489, (49, 2) = -29.39558090446206, (50, 1) = 11.79378728978802, (50, 2) = -40.669372493045, (51, 1) = 11.183836634124887, (51, 2) = -48.03607321421275, (52, 1) = 10.502833198232775, (52, 2) = -51.72939943874306, (53, 1) = 9.796155822525524, (53, 2) = -52.452900134134595, (54, 1) = 9.096522907410908, (54, 2) = -51.04757024201683, (55, 1) = 8.427075143053173, (55, 2) = -48.28611154332724, (56, 1) = 7.801521511288004, (56, 2) = -44.78604220693657, (57, 1) = 7.22588851589702, (57, 2) = -40.986432776247455, (58, 1) = 6.701837954309735, (58, 2) = -37.176690679157566, (59, 1) = 6.228168238234758, (59, 2) = -33.53206368300497}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(59, {(1) = .0, (2) = 0.13680812436809657e-1, (3) = 0.27971051145644167e-1, (4) = 0.4296023135974479e-1, (5) = 0.5851786182852681e-1, (6) = 0.7435973609701907e-1, (7) = 0.9049963296106847e-1, (8) = .1070194360291177, (9) = .12400404122395811, (10) = .14149911102616744, (11) = .15962080170244286, (12) = .17846661163773977, (13) = .1981659899247577, (14) = .2189484735830971, (15) = .2410342106109193, (16) = .2648591917099649, (17) = .29096256417430655, (18) = .3201405168424954, (19) = .3537506263213988, (20) = .39359925436482457, (21) = .44328213017143814, (22) = .508829664760916, (23) = .5993459008179163, (24) = .7156934169226271, (25) = .8282966673508598, (26) = .9130409367878883, (27) = .9748398788122531, (28) = 1.0223008015052495, (29) = 1.0606661846616776, (30) = 1.0932188074676115, (31) = 1.1216707736458404, (32) = 1.1472399927159063, (33) = 1.1706725764522086, (34) = 1.1924372702961963, (35) = 1.2129338175389266, (36) = 1.2323858241405192, (37) = 1.250955152949817, (38) = 1.268795762648431, (39) = 1.2859897085262963, (40) = 1.3026146566339305, (41) = 1.3187452976919192, (42) = 1.3344237627757964, (43) = 1.3497030677053556, (44) = 1.3646567912885574, (45) = 1.3793073520205892, (46) = 1.3936816310415334, (47) = 1.407863413800478, (48) = 1.421874455055364, (49) = 1.4357209480594326, (50) = 1.449452395680341, (51) = 1.463105461314161, (52) = 1.4766814720936998, (53) = 1.4901926102416252, (54) = 1.5036746366401115, (55) = 1.5171303716739315, (56) = 1.5305608572874234, (57) = 1.5439797161303026, (58) = 1.5573915356105725, (59) = 1.5707963267949}, datatype = float[8], order = C_order); Y := Matrix(59, 2, {(1, 1) = .0, (1, 2) = -0.2040971830865995e-12, (2, 1) = -0.2515085127931946e-15, (2, 2) = -0.2182240868305703e-12, (3, 1) = -0.3034400215680231e-14, (3, 2) = -0.24629572182126297e-12, (4, 1) = -0.70293294530316345e-14, (4, 2) = -0.27355561891671303e-12, (5, 1) = -0.11642509457266535e-13, (5, 2) = -0.2921725988376451e-12, (6, 1) = -0.14891698174826332e-13, (6, 2) = -0.27332966679506926e-12, (7, 1) = -0.2244158402930937e-13, (7, 2) = -0.21555448086863683e-12, (8, 1) = -0.23830070405255755e-13, (8, 2) = -0.14452163468605652e-12, (9, 1) = -0.28784241689197694e-13, (9, 2) = -0.8936429016769738e-13, (10, 1) = -0.2870763992607896e-13, (10, 2) = -0.7015153701410498e-13, (11, 1) = -0.29805232120338505e-13, (11, 2) = -0.20455226198960395e-13, (12, 1) = -0.3071869736898188e-13, (12, 2) = -0.949608712395301e-14, (13, 1) = -0.33085750279284234e-13, (13, 2) = -0.9853725522439467e-14, (14, 1) = -0.3024025650769599e-13, (14, 2) = 0.55343272382320285e-14, (15, 1) = -0.2989762453456761e-13, (15, 2) = 0.9484091577415516e-14, (16, 1) = -0.3126420836975352e-13, (16, 2) = 0.10993914919372806e-15, (17, 1) = -0.2904952409517121e-13, (17, 2) = 0.7741875869166384e-14, (18, 1) = -0.26389405390354672e-13, (18, 2) = -0.39736138956257974e-14, (19, 1) = -0.3066642976676729e-13, (19, 2) = -0.8841791351636723e-14, (20, 1) = -0.3381874764737987e-13, (20, 2) = -0.10474029579197825e-13, (21, 1) = -0.3076950762715862e-13, (21, 2) = -0.14630357321113296e-13, (22, 1) = -0.25742990878046857e-13, (22, 2) = -0.32013707096004785e-13, (23, 1) = -0.26034531106847924e-14, (23, 2) = -0.17882398235769704e-13, (24, 1) = -0.24294629294277565e-13, (24, 2) = 0.6033751098446911e-13, (25, 1) = -0.10212795741481291e-13, (25, 2) = -0.5409909790245115e-14, (26, 1) = 0.65234775286404595e-14, (26, 2) = 0.10506976038853487e-14, (27, 1) = 0.15145921627767115e-13, (27, 2) = 0.11062582034309036e-13, (28, 1) = 0.17727034515154246e-13, (28, 2) = 0.1392374003706416e-13, (29, 1) = 0.20190304444405087e-13, (29, 2) = 0.24077662756883795e-13, (30, 1) = 0.15567756982704383e-13, (30, 2) = 0.3220356728033926e-13, (31, 1) = 0.2049236289955625e-13, (31, 2) = 0.4582315826902611e-13, (32, 1) = 0.2093779309445625e-13, (32, 2) = 0.50847503530540665e-13, (33, 1) = 0.34656722670579105e-13, (33, 2) = 0.65791408284367e-13, (34, 1) = 0.1414627403891834e-13, (34, 2) = 0.7434970442839986e-13, (35, 1) = 0.1430197385087178e-13, (35, 2) = 0.807995248637132e-13, (36, 1) = 0.2674246752632126e-13, (36, 2) = 0.8238316242997303e-13, (37, 1) = 0.3749521089102351e-13, (37, 2) = 0.9551972331079688e-13, (38, 1) = 0.43798173157104124e-13, (38, 2) = 0.8882243056529444e-13, (39, 1) = 0.34765247192333954e-13, (39, 2) = 0.7174268116999525e-13, (40, 1) = 0.3492718539460031e-13, (40, 2) = 0.25074509668472094e-13, (41, 1) = 0.3268343950774218e-13, (41, 2) = -0.6182099171300564e-13, (42, 1) = 0.7615127467050659e-14, (42, 2) = -0.1854246486047524e-12, (43, 1) = 0.20280258788197414e-13, (43, 2) = -0.2612655964015281e-12, (44, 1) = 0.5268977580072408e-14, (44, 2) = -0.3058713868413555e-12, (45, 1) = 0.8243135377351008e-14, (45, 2) = -0.3000324869932838e-12, (46, 1) = 0.9651203882669847e-14, (46, 2) = -0.23958986839330786e-12, (47, 1) = -0.15172603499122148e-13, (47, 2) = -0.14114952424844551e-12, (48, 1) = 0.6415293117671377e-14, (48, 2) = -0.3574626525598719e-13, (49, 1) = 0.11771473310921094e-13, (49, 2) = 0.12758808412240534e-13, (50, 1) = 0.350361321759921e-14, (50, 2) = 0.18960251337597162e-13, (51, 1) = -0.402567915071639e-14, (51, 2) = -0.42598728779319296e-13, (52, 1) = 0.14493706543893737e-13, (52, 2) = -0.4380811741884113e-13, (53, 1) = 0.19531436506149837e-13, (53, 2) = -0.8490114617635297e-13, (54, 1) = 0.4784180312117899e-14, (54, 2) = -0.10238102255674353e-12, (55, 1) = 0.13186433513432542e-13, (55, 2) = -0.12142340951879697e-12, (56, 1) = 0.68589188913736246e-14, (56, 2) = -0.10734720880680753e-12, (57, 1) = 0.4428481028761533e-14, (57, 2) = -0.9866461770792047e-13, (58, 1) = 0.7387414481358141e-14, (58, 2) = -0.9166293971010342e-13, (59, 1) = .0, (59, 2) = -0.731647525953621e-13}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[59] elif outpoint = "order" then return 10 elif outpoint = "error" then return HFloat(3.058713868413555e-13) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 59, [v(u), diff(v(u), u)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[59] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[59] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(59, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(59, 2, X, Y, outpoint, yout, L, V) end if; [u = outpoint, seq('[v(u), diff(v(u), u)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[59] elif outpoint = "order" then return 10 elif outpoint = "error" then return HFloat(3.058713868413555e-13) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 59, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[59] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[59] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(59, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(59, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(1..3, {(1) = 36893628625353841228, (2) = 36893628625353841404, (3) = 36893628625353841580}), (3) = [u, v(u), diff(v(u), u)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(u) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(u) else `D(v)` := pointto(data[2][3]); return ('`D(v)`')(u) end if end if; try res := solnproc(outpoint); res[3] catch: error  end try end proc]

display(
        torus,
        tubeplot([seq(eval(x(u, v(u)), dsol1))], u=0..Pi/2, color=red, radius=0.1),
        tubeplot([seq(eval(x(u, v(u)), dsol3))], u=0..Pi/2, color=cyan, radius=0.1),
        pointplot3d(eval([x(0,v(0)),x(Pi/2,v(Pi/2))], dsol3), symbol=solidsphere, symbolsize=20, color="Orange"),
lightmodel=light4, thickness=5, orientation=[-65,45,0], size=[700,400],
style=surface, axes=none, viewpoint=circleright);

 

 

>

 

Download geodesics-between-two-points.mw

You will find the geodesic equations in the wikipedia web page that you have referenced to.   Admittedly it is rather cryptic, involving Christoffel symbols in an implicit summation.  Once we unpack the notation, the geodesic equations for a surface turn out to be a pair of coupled, second order, nonlinear ODEs.  These equations are usually derived in an introductory course on differential geometry.  To get to those equations, one needs to learn about a surface's first fundamental form and how the Christoffel symbols are expressed in terms of the coefficients of the first fundamental form.  

The Christoffel() proc in the attached worksheet calculates the coefficients of the first fundamental form of any surface and then calculates and returns the corresponding Christoffel symbols.  To understand the details, you need to know some differential geometry.

A second proc, Geodesic_DEs(), applies the equation in the Wikipedia page to calculate and return the surface's geodesic equations.  Once you have those equations, you can apply Maple's dsolve() to calculate and plot the geodesics.  The worksheet includes an illustrative example.

Maple's Differential Geometry package has a GeodesicEquations() function.  Odds are that what I have calculated here can be done more quickly by applying that function but don't know how.  Perhaps someone who knows may show us.

Answers to the several questions that you asked would involve some elementary differential geometry, so I cannot summarize them here in a few lines.  I can, however, give short answers to two of those questions:

  1. What does "to walk forward" mean?  When you move on a surface, at each location you have a velocity vector, an acceleration vector, and a normal vector (i.e., vector perpendicular to the surface).  You are "moving forward", or equivalently, "traveling on a geodesic", provided that those three vectors are co-planar.
  2. Geodesic is the shortest distance between two points in space, but no such end points are shown in the website's animation.  A geodesic typically extends indefinitely in both directions.  There are no end points.  A curve is geodesic if for any pair of points A, B on the curve, the part of the curve that lies between A and B is the shortest path among all curves that lie on the surface and connect A to B. 

restart;

with(plots):

Receives a parametrization x of a surface and returns the corresponding

Christoffel symbols `#msubsup(mi("&Gamma;",fontstyle = "normal"),mi("jk"),mi("i"))`, each of which is a procedure.

Christoffel := proc(x::procedure)
        local u, v, xu, xv, E, F, G, Eu, Ev, Fu, Fv, Gu, Gv, den, C;
        xu := diff(x(u,v),u);
        xv := diff(x(u,v),v);
        # E, F, G are the coefficients of the surface's first fundamental form
        E := xu^+ . xu;
        F := xu^+ . xv;
        G := xv^+ . xv;
        Eu := diff(E, u);
        Ev := diff(E, v);
        Fu := diff(F, u);
        Fv := diff(F, v);
        Gu := diff(G, u);
        Gv := diff(G, v);
        den := 2*(E*G - F^2);
        C[1][1,1] := unapply(simplify((G*Eu - 2*F*Fu + F*Ev)/den), u, v);
        C[1][1,2] := unapply(simplify((G*Ev - F*Gu)/den), u, v);
        C[1][2,2] := unapply(simplify(-(F*Gv - 2*G*Fv + G*Gu)/den), u, v);
        C[2][1,1] := unapply(simplify(-(E*Ev - 2*E*Fu + F*Eu)/den), u, v);
        C[2][1,2] := unapply(simplify((E*Gu - F*Ev)/den), u, v);
        C[2][2,2] := unapply(simplify((E*Gv - 2*F*Fv + F*Gu)/den), u, v);
        C[1][2,1] := C[1][1,2];
        C[2][2,1] := C[2][1,2];
        return C;
end proc:

Example: A sphere

x := (u,v) -> < cos(v)*cos(u), cos(v)*sin(u), sin(v) >;

proc (u, v) options operator, arrow; `<,>`(cos(v)*cos(u), cos(v)*sin(u), sin(v)) end proc

Here are the 8 Christoffel symbols of our parametrization of the sphere:

Gamma := Christoffel(x):

Gamma[1][1,1](u,v), Gamma[1][1,2](u,v), Gamma[1][2,1](u,v), Gamma[1][2,2](u,v);

0, -sin(v)/cos(v), -sin(v)/cos(v), 0

Gamma[2][1,1](u,v), Gamma[2][1,2](u,v), Gamma[2][2,1](u,v), Gamma[2][2,2](u,v);

cos(v)*sin(v), 0, 0, 0

Example: A torus of major radius a and minor radius bNULL

x := (u,v) -> < (a+b*cos(v))*cos(u), (a+b*cos(v))*sin(u),  b*sin(v) >;

proc (u, v) options operator, arrow; `<,>`((a+b*cos(v))*cos(u), (a+b*cos(v))*sin(u), b*sin(v)) end proc

Gamma := Christoffel(x):

Gamma[1][1,1](u,v), Gamma[1][1,2](u,v), Gamma[1][2,1](u,v), Gamma[1][2,2](u,v);

0, -b*sin(v)/(a+b*cos(v)), -b*sin(v)/(a+b*cos(v)), 0

Gamma[2][1,1](u,v), Gamma[2][1,2](u,v), Gamma[2][2,1](u,v), Gamma[2][2,2](u,v);

(a+b*cos(v))*sin(v)/b, 0, 0, 0

The geodesic equations

Geodesic_DEs := proc(x::procedure)
        local Gamma, de, i, j, k, X;
        Gamma := Christoffel(x);
        for i from 1 to 2 do
           de[i] := diff(X[i](t), t, t) + add(add(
                        Gamma[i][j,k](X[1](t),X[2](t))*diff(X[j](t),t)*diff(X[k](t),t),
                        j=1..2), k=1..2) = 0;
        end do:
        eval([de[1],de[2]], {X[1]=u, X[2]=v})[];   # change notation from X to u, v
        return (%);
end proc:

 

Geodesics on a torus

 

Here we calculate and plot geodesics on a torus.  You may change the torus

to any other parametrized surface.  Everthing else remains the same.

< (a+b*cos(v))*cos(u), (a+b*cos(v))*sin(u),  b*sin(v) >;
eval(%, {a=5, b=2});
x := unapply(%, u, v):

Vector(3, {(1) = (a+b*cos(v))*cos(u), (2) = (a+b*cos(v))*sin(u), (3) = b*sin(v)})

Vector[column](%id = 36893628224392653508)

the_torus := plot3d(x(u,v), u=-Pi..Pi, v=-Pi..Pi, scaling=constrained, style=wireframe);

The differential equations of geodesics on torus

DEs := Geodesic_DEs(x);

diff(diff(u(t), t), t)-4*sin(v(t))*(diff(v(t), t))*(diff(u(t), t))/(5+2*cos(v(t))) = 0, diff(diff(v(t), t), t)+(cos(v(t))*sin(v(t))+(5/2)*sin(v(t)))*(diff(u(t), t))^2 = 0

Pick any desired starting point and direction on the torus.  The angle theta determines
the geodesic's initial direction.

unassign('theta');
ICs := u(0) = 0, v(0) = 0, D(u)(0) = cos(theta), D(v)(0) = sin(theta);

u(0) = 0, v(0) = 0, (D(u))(0) = cos(theta), (D(v))(0) = sin(theta)

Pick a theta and solve the geodesic DEs

theta:= 0.95*Pi/2:
dsol := dsolve({DEs, ICs}, numeric, output=operator);