Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

The try/catch mechanism may be what you are looking for.

restart;

foo:=proc(n::integer)
        print("before internal proc");
        proc()
                if n<0 then
                        error "something";
                fi;
        end proc(); #notice the (); at end
        print("after internal proc");
end proc:

Let's try:

try foo(3); (a+b)^2;  catch "something": (c+d)^2; end try;

"before internal proc"

"after internal proc"

(a+b)^2

try foo(-3); (a+b)^2;  catch "something": (c+d)^2; end try;

"before internal proc"

(c+d)^2

 

If you feel very enthusiastic, you may bundle up the try/catch along with foo
into a signle proc.

restart;

bar := proc(n::integer)
        local foo;
        foo:=proc(n::integer)
                print("before internal proc");
                proc()
                        if n<0 then
                                error "something";
                        fi;
                end proc(); #notice the (); at end
                print("after internal proc");
        end proc:
        try foo(n); (a+b)^2;  catch "something": (c+d)^2; end try;
end proc:

 

Let's try:

bar(3);

"before internal proc"

"after internal proc"

(a+b)^2

bar(-3);

"before internal proc"

(c+d)^2

Download mw.mw

Your code expresses the essence of the solution, but it can be streamlined a bit.

restart;
eq := <x, y> =~ lambda*<1, 1> + mu*<1, 2>;
convert(eq, set);
solve(%, {lambda, mu});

We see that solve() returns a unique solution, indicating that any vector <x,y> may be uniquely represented as a linear combination of <1,1> and <1,2>, indicating that those form a basis for the space.

Maple's PDF export has been broken in the last several releases. Despite repeated complaints on this forum, no solution seems to be forthcoming.

These two workarounds work quite well on the command-line in Linux.  The equivalent command lines may be available in Mac and Windows but I have neither so I wouldn't know.

  1. Workaround #1:  Save the graphics as EPS, let's say filename.eps, then do
        epstopdf filename.eps
    to produce filename.pdf.
  2. Workaround #2:  Save the graphics as SVG, let's say filename.svg, then do
        inkscape --export-type=pdf --export-area-drawing filename.svg
    to produce filename.pdf.

Both methods produce good quality PDF.  The file size of the PDF file in method #2 is about 1/20 of that of method #1 but the latter is more faithful to the original.

I have attached PDFs cooresponding to 
    plot3d(x^2+y^2, x=-1..1, y=-1..1);
converted through the two methods.

parabola-from-eps.pdf
parabola-from-svg.pdf

 

Here is the solution, but don't just pass it on to your teacher. She won't believe you.

To get an answer with detailed explanation, show what you know and how far you have gotten on your own on this homework.

restart;
f := x -> (x^3 + 9*x^2 - 9*x - 1) / (x^4 + 1);
numer(D(f)(x) - 1/2);
xvals := fsolve(%, x, maxsols=4);
seq(f(t) + 0.5*(x-t), t in xvals);
plot([f(x), %], x=-7..7, scaling=constrained, color=[seq(cat("oldplots ", i), i=1..5)]);

restart;
with(plots):
with(plottools):
p1 := plot(sin(x), x=-2*Pi..2*Pi, color="Red"):
display(
	plot(sin(x), x=-2*Pi..2*Pi, color="Green"),
	line([-2*Pi,0], [2*Pi,0]),
	line([0,-2*Pi], [0,2*Pi])
):
p2 := rotate(%, Pi/4):
display(p1,p2, scaling=constrained);

I don't see a simple way of adding tickmarks.

restart;
with(plots):
vertices := [
	c__1 = <0, -2>,
	c__2 = <1, 2>,
	c__3 = <2, 2>,
	c__4 = <0.5, 6>,
	c__5 = <1, 4> 
]:
Labels := map( v -> [convert(rhs(v), list)[], lhs(v)], vertices):
display(
	textplot(Labels, align={above}),
	plottools:-polygon(convert~(rhs~(vertices), list), style=line, color=red)
);

 

See https://www.mapleprimes.com/questions/229945-Simulation-Of-A-Solution-Of--Neutral for a lot of ideas on how to handle delay differential equations.

 

I don't see what you mean by tracing the moving points since the number of point changes in your sequence of plots. Here is how to draw traces of individual moving points that preserve their identities during the motion.

restart;
with(plots):
f := t -> sin(t);
g := t -> [cos(t), sin(t)];
trace_me := proc(t)
	display(
		pointplot([t,f(t)], symbol=solidcircle, symbolsize=30, color="Red"),
		plot(f(s), s=0..t, color="Green"),
		pointplot(g(t), symbol=solidcircle, symbolsize=30, color="Orange"),
		plot([g(s)[], s=0..t], color="Blue"),
	scaling=constrained);
end proc:
animate(trace_me, [t], t=0..2*Pi);

I will illustrate the required procedure in terms of a function that I made up. If this does not work in your case, then you need to provide more details about your function fp(x).

restart;

with(plots):

Here is a family of function of x defined in terms of a parameter p.

f := (x,p) -> (x-p)^2 + cos(p);

proc (x, p) options operator, arrow; (x-p)^2+cos(p) end proc

To calculate the family's envelope, we solve the following system
of two equations for the two unknowns x and y.  Call the solution X(p)

and Y(p).  The envelope is the parametric curve X(p), Y(p).

y = f(x,p), diff(f(x,p),p) = 0;
solve({%}, {x,y}):
X := unapply(eval(x, %), p);
Y := unapply(eval(y, %%), p);

y = (x-p)^2+cos(p), -2*x+2*p-sin(p) = 0

proc (p) options operator, arrow; p-(1/2)*sin(p) end proc

proc (p) options operator, arrow; (1/4)*sin(p)^2+cos(p) end proc

Let us plot the result.  The envelope is shown in blue.

display(
        seq(plot(f(x,p), x=-10..10, color="Green", view=-1..2), p=-8..8, 0.3),
        plot([X(p),Y(p), p=-8..8], color="Blue", thickness=3),
size=[700,500]);

We now calculate the area bounded by the central arch of the envelope and the x axis.
From the graph above we see that the arch extends roughly from x = -1 to x = 1.

To calculate the precise range, we seek parameter values that make Y(p) = 0

p_roots := solve(Y(p), p);

arctan(2*(-2+5^(1/2))^(1/2)/(-5^(1/2)+2))+Pi, -arctan(2*(-2+5^(1/2))^(1/2)/(-5^(1/2)+2))-Pi, arctan(2*(-2-5^(1/2))^(1/2), 5^(1/2)+2), arctan(-2*(-2-5^(1/2))^(1/2), 5^(1/2)+2)

and then calculate the corresponding x values:

simplify(X~([p_roots]));
evalf(%);

[-arctan(2/(-2+5^(1/2))^(1/2))+Pi-(-2+5^(1/2))^(1/2), arctan(2/(-2+5^(1/2))^(1/2))-Pi+(-2+5^(1/2))^(1/2), arctan((2*I)*(5^(1/2)+2)^(1/2), 5^(1/2)+2)-I*(5^(1/2)+2)^(1/2), arctan(-(2*I)*(5^(1/2)+2)^(1/2), 5^(1/2)+2)+I*(5^(1/2)+2)^(1/2)]

[1.323245518, -1.323245518, 0.64379097e-1*I, -0.64379097e-1*I]

We conclude that the support of the arch is roughly in the range abs(x) < 1.323. The corresponding

parameter values are

p_roots[1], p_roots[2];

arctan(2*(-2+5^(1/2))^(1/2)/(-5^(1/2)+2))+Pi, -arctan(2*(-2+5^(1/2))^(1/2)/(-5^(1/2)+2))-Pi

The area bounded by a parametric curve and the x axis is given by int(Y(p)*(D(X))(p), p):

A := simplify(int(Y(p)*D(X)(p), p = p_roots[2] .. p_roots[1]));

(7/3)*(-2+5^(1/2))^(1/2)+(5/6)*(-2+5^(1/2))^(1/2)*5^(1/2)+(1/4)*arctan(2*(-2+5^(1/2))^(1/2)*(5^(1/2)+2))-(1/4)*Pi

evalf(A);

1.586776254
 

Download area-under-envelope.mw

You are interpreting Maple's result as if it had said y(x) = whatever.  But that's not what the answer says.  Rather, the answer is an integral with the upper limit as y(x).  The following worksheet shows how to verify that Maple's solution does indeed satisfy the differential equation.

restart;

ode := diff(y(x), x, x) = - sin(y(x));   # the minus sign is missing in the original worksheet!

diff(diff(y(x), x), x) = -sin(y(x))

dsol := dsolve(ode);

Intat(1/(2*cos(_a)+_C1)^(1/2), _a = y(x))-x-_C2 = 0, Intat(-1/(2*cos(_a)+_C1)^(1/2), _a = y(x))-x-_C2 = 0

Let us verify that dsol[1] is a solution of the ode.  (dsol[2] may be verified in the same way)

diff(dsol[1], x);
isolate(%, diff(y(x), x));
diff(%, x);
applyop(eval, 2, %, %%);

(diff(y(x), x))/(2*cos(y(x))+_C1)^(1/2)-1 = 0

diff(y(x), x) = (2*cos(y(x))+_C1)^(1/2)

diff(diff(y(x), x), x) = -(diff(y(x), x))*sin(y(x))/(2*cos(y(x))+_C1)^(1/2)

diff(diff(y(x), x), x) = -sin(y(x))

 

 

Download mw.mw

If we plot the left-hand side and right-hand side expressions of your equation, we see that the plots don't intersect, and therefore there is no (real) solution to your equation.

eq := exp(x) * (1 + sqrt(1 + exp(-2*x))) - arcsinh(exp(-x)) = x;
plot([lhs(eq), rhs(eq)], x=-3..3);

Here is how to go about it.  The resulting figure is not exactly like what you have shown which indicates that your figure was produced based on a different set of parameters than what you have provided.

restart;

We have five differential equations:

de1 := diff(u__1(y),y,y) + 1/2*diff(N(y),y) + 1/2*theta__1(y) = 0;

diff(diff(u__1(y), y), y)+(1/2)*(diff(N(y), y))+(1/2)*theta__1(y) = 0

de2 := diff(N(y),y,y) - 2/3 * (2*N(y) + diff(u__1(y),y)) = 0;

diff(diff(N(y), y), y)-(4/3)*N(y)-(2/3)*(diff(u__1(y), y)) = 0

de3 := diff(theta__1(y),y,y) = 0;

diff(diff(theta__1(y), y), y) = 0

de4 := diff(u__2(y),y,y) + theta__2(y) = 0;

diff(diff(u__2(y), y), y)+theta__2(y) = 0

de5 := diff(theta__2(y),y,y) = 0;

diff(diff(theta__2(y), y), y) = 0

Find the general solution:

dsol := dsolve({de1,de2,de3,de4,de5});

{N(y) = exp(y)*_C2-exp(-y)*_C3+(1/3)*_C10*y+(1/6)*_C9*y^2+_C4, u__1(y) = -2*_C4*y-(1/3)*_C10*y^2-(1/9)*_C9*y^3-(1/2)*exp(y)*_C2-(1/2)*exp(-y)*_C3+(1/2)*_C9*y+_C1, u__2(y) = -(1/6)*_C7*y^3-(1/2)*_C8*y^2+_C5*y+_C6, theta__1(y) = _C9*y+_C10, theta__2(y) = _C7*y+_C8}

How many integrations constants?

indets(dsol, name);

{_C1, _C10, _C2, _C3, _C4, _C5, _C6, _C7, _C8, _C9, y}

So there are 10 integration constants, which is good since we have 10 boundary and interface conditions.

 

For now, let us extract the individual solutions of the ODE system:

u__1 := 'u__1':
eval(u__1(y), dsol):
u__1 := unapply(%, y);

proc (y) options operator, arrow; -2*_C4*y-(1/3)*_C10*y^2-(1/9)*_C9*y^3-(1/2)*exp(y)*_C2-(1/2)*exp(-y)*_C3+(1/2)*_C9*y+_C1 end proc

u__2 := 'u__2':
eval(u__2(y), dsol):
u__2 := unapply(%, y);

proc (y) options operator, arrow; -(1/6)*_C7*y^3-(1/2)*_C8*y^2+_C5*y+_C6 end proc

theta__1 := 'theta__1':
eval(theta__1(y), dsol):
theta__1 := unapply(%, y);

proc (y) options operator, arrow; _C9*y+_C10 end proc

theta__2 := 'theta__2':
eval(theta__2(y), dsol):
theta__2 := unapply(%, y);

proc (y) options operator, arrow; _C7*y+_C8 end proc

N := 'N':
eval(N(y), dsol):
N := unapply(%, y);

proc (y) options operator, arrow; exp(y)*_C2-exp(-y)*_C3+(1/3)*_C10*y+(1/6)*_C9*y^2+_C4 end proc

Now let us look at the boundary condtions.

 

Boundary conditions on the interval (-1,0):

bc[-1] := u__1(-1) = 0, N(-1) = 0, theta__1(-1) = 1,
        u__1(0) = u__0, theta__1(0) = theta__0, D(theta__1)(0) = theta_d;

2*_C4-(1/3)*_C10-(7/18)*_C9-(1/2)*exp(-1)*_C2-(1/2)*exp(1)*_C3+_C1 = 0, exp(-1)*_C2-exp(1)*_C3-(1/3)*_C10+(1/6)*_C9+_C4 = 0, -_C9+_C10 = 1, -(1/2)*_C2-(1/2)*_C3+_C1 = u__0, _C10 = theta__0, _C9 = theta_d

Boundary conditions on the interval (0,1)

bc[1] := u__2(1) = 0, theta__2(1) = 0;

-(1/6)*_C7-(1/2)*_C8+_C5+_C6 = 0, _C7+_C8 = 0

Interface conditions at  y = 0

bc[0] := u__1(0) = u__2(0), D(N)(0) = 0, theta__1(0) = theta__2(0),  D(theta__1)(0) = D(theta__2)(0),
        D(u__1)(0) + 1/2*N(0) = 1/2*D(u__2)(0);

-(1/2)*_C2-(1/2)*_C3+_C1 = _C6, _C2+_C3+(1/3)*_C10 = 0, _C10 = _C8, _C9 = _C7, -(3/2)*_C4+(1/2)*_C9 = (1/2)*_C5

Apply the 10 boundary and interface conditions to determine the 10 integration coefficients:

c_vals := solve({bc[-1], bc[0], bc[1]});

{_C1 = -(1/6)*(6*exp(-1)*exp(1)-21*exp(-1)-5*exp(1))/(11*exp(-1)+9*exp(1)), _C10 = 1/2, _C2 = -(1/18)*(27*exp(1)-67)/(11*exp(-1)+9*exp(1)), _C3 = -(1/18)*(33*exp(-1)+67)/(11*exp(-1)+9*exp(1)), _C4 = -(1/36)*(12*exp(-1)*exp(1)+35*exp(-1)+53*exp(1))/(11*exp(-1)+9*exp(1)), _C5 = (1/12)*(12*exp(-1)*exp(1)-31*exp(-1)-exp(1))/(11*exp(-1)+9*exp(1)), _C6 = -(1/12)*(12*exp(-1)*exp(1)-53*exp(-1)-19*exp(1))/(11*exp(-1)+9*exp(1)), _C7 = -1/2, _C8 = 1/2, _C9 = -1/2, theta_d = -1/2, u__0 = -(1/12)*(12*exp(-1)*exp(1)-53*exp(-1)-19*exp(1))/(11*exp(-1)+9*exp(1)), theta__0 = 1/2}

Evaluate u1 and u2

my_u1 := eval(u__1(y), c_vals);
my_u2 := eval(u__2(y), c_vals);

(1/18)*(12*exp(-1)*exp(1)+35*exp(-1)+53*exp(1))*y/(11*exp(-1)+9*exp(1))-(1/6)*y^2+(1/18)*y^3+(1/36)*exp(y)*(27*exp(1)-67)/(11*exp(-1)+9*exp(1))+(1/36)*exp(-y)*(33*exp(-1)+67)/(11*exp(-1)+9*exp(1))-(1/4)*y-(1/6)*(6*exp(-1)*exp(1)-21*exp(-1)-5*exp(1))/(11*exp(-1)+9*exp(1))

(1/12)*y^3-(1/4)*y^2+(1/12)*(12*exp(-1)*exp(1)-31*exp(-1)-exp(1))*y/(11*exp(-1)+9*exp(1))-(1/12)*(12*exp(-1)*exp(1)-53*exp(-1)-19*exp(1))/(11*exp(-1)+9*exp(1))

plots:-display(
        plot(my_u1, y=-1..0, color=red),
        plot(my_u2, y=0..1, color=blue)
);

 

Download mw.mw

 

As I noted earlier, your equation is linear, so your focus on nonlinear equations is out of place. If you refer to the discussion of linear first order equations in your differential equations textbook, you will find a general solution obtained through the method of integrating factors. The transient and steady-state parts of the solution may be extracted from there.

What you should be asking is: What do the transient and steady-state solutions look like? I have attached a PDF in which I have analyzed a slightly more general differential equation than the one that you have asked about. In the following worksheet I use the results obtained in that PDF to calculate and plot the transient and steady-state graphs for your equation. Note that this has nothing to do with nonlinear equations.

restart;

with(plots):

In the attached PDF I have analyzed the following first order linear equation:

diff(x(t),t) = p(t) - (a + q(t))*x(t);

diff(x(t), t) = p(t)-(a+q(t))*x(t)

Assumption:

• 

p(t) and q(t) are bounded periodic functions of some period, say σ.

• 

the average of q(t) over a period is zero

• 

a > 0 is a constant.

Define

Q(t) = Int(q(s), s=0..t);

Q(t) = Int(q(s), s = 0 .. t)

It is evident that Q is σ-periodic and Q(0) = 0.


In the attached PDF I show that the solution x(t) subject to an initial condition x(0)=x_0

tends to a σ-periodic steady state as time goes to infinity.  The transient and steady state

parts of the solution are given by:

x__tr(t) = x__0*exp(-a*t-Q(t));
x__ss(t) = Int(p(t-s)*exp(-a*s - Q(t) + Q(t-s)), s=0..T);

x__tr(t) = x__0*exp(-a*t-Q(t))

x__ss(t) = Int(p(t-s)*exp(-a*s-Q(t)+Q(t-s)), s = 0 .. T)

where the constant T is a sufficiently large time value beyond which

the transient may be considered as died out.

 

A concrete example

 

de := diff(x(t),t) = 1 - ( 1 + sin(10*t)/10) * x(t), x(0) = 0.9;

diff(x(t), t) = 1-(1+(1/10)*sin(10*t))*x(t), x(0) = .9

This corresponds to the choices of p(t) = 1, q(t) = sin(10*t)/10, and a = 1. We calculate Q:

int(sin(10*s)/10, s=0..t):
Q := unapply(%, t);

proc (t) options operator, arrow; 1/100-(1/100)*cos(10*t) end proc

Let us solve the differential equation numerically by applying Maple's dsolve():

dsol := dsolve({de}, numeric);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 27, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..65, {(1) = 1, (2) = 1, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0, (64) = -1, (65) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = .5047658755841546, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..1, {(1) = .9}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..1, {(1) = .9}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0.9999999999999998e-1}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; YP[1] := 1-(1+(1/10)*sin(10*X))*Y[1]; 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = x(t)]`; YP[1] := 1-(1+(1/10)*sin(10*X))*Y[1]; 0 end proc, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 27 ) = (""), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..1, {(1) = 0.}); _vmap := array( 1 .. 1, [( 1 ) = (1)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch:  end try else try procname(procname("right")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch:  end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error  end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, x(t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error  end try end proc

and graph it:

p__exact := odeplot(dsol, t=0..10, color=red);

We calculate the transient from equation (4a) in the attached PDF and plot it:

x__tr := 0.9*exp(-t-Q(t));

.9*exp(-t-1/100+(1/100)*cos(10*t))

plot(x__tr, t=0..10);

We see that the transient dies out at about T=7, beyond which the solution

goes into steady-state. The steady-state is given in equation (4b) in the PDF:

x__ss := Int(exp(-s-Q(t) + Q(t-s)), s=0..7);

Int(exp(-s+(1/100)*cos(10*t)-(1/100)*cos(-10*t+10*s)), s = 0 .. 7)

and here is what it looks like:

p__ss := plot(x__ss, t=0..10, color=blue);

Let's plot the exact solution and the steady-state together to verify

that they match:

display(p__exact, p__ss);

 

 

 

Download worksheet: steady-state-of-an-ode-2.mw

Download documentation: steady-state-of-an-ode.pdf

Your equation has two solutions, and therefore sol produced by maple consists of two parts:

sol := {solution1}, {solution2}

In your code, you are attempting to convert sol to Matlab() and that backfires.  You need to pass the individual solution1 and solution2 to Matlab():

Matlab(sol[1][]);
Matlab(sol[2][]);

The purpose of the empty brackets is to remove the curly braces in sol.

 

The "triangle" that you have depicted is not a spherical triangle since the blue edge is not an arc of a great circle, and therefore the Pythagorean Theorem does not apply.

Aside: You have not shown your calculations, so I calculated everything from scratch.  My calculation shows that the length of the green curve is 

which evaluates to 0.9228460424.

1 2 3 4 5 6 7 Last Page 2 of 44