Advanced Search

Maple Questions and Posts

These are Posts and Questions associated with the product, Maple
View unanswered posts
Maple Questions and Posts Feed

is(0 < x^2 + x + y) assuming (0 < x + y) wrongly g...

Asked by: tilo rc 25
assumptions assuming
February 16 2025
1  1

"is(0 < x^2 + x + y) assuming (0 < x + y)" wrongly returns false. I'm using Maple 2024.2. This is a bug, right?

How collect a fraction which contain derivative? A...

Asked by: salim-barzani 655
parameters pde differential-equations
February 15 2025
2  10

i want construct a series trail function for all pdf not just this one but this is a easy one, also after replacing the function How i can collect variable and make algebraic system for finding the constant of series function like a[20],a[10],a[00]. where i is number of derivative by x and n is number of derivative by t also n=0 and m=2 

> 

restart
> 

with(PDEtools)
> 

with(LinearAlgebra)
> 

NULL
> 

with(SolveTools)
> 

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

 (1)
> 

declare(u(x, t))

u(x, t)*`will now be displayed as`*u

 (2)
> 

declare(w(x, t))

w(x, t)*`will now be displayed as`*w

 (3)
> 

pde := diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x))+delta*(diff(u(x, t), `$`(x, 3))) = 0

diff(u(x, t), t)+u(x, t)*(diff(u(x, t), x))+delta*(diff(diff(diff(u(x, t), x), x), x)) = 0

 (4)
> 

NULL
> 

K := u(x, t) = a[20]*(diff(ln(w(x, t)), `$`(x, 2)))+a[10]*(diff(ln(w(x, t)), x))+a[0]

u(x, t) = a[20]*((diff(diff(w(x, t), x), x))/w(x, t)-(diff(w(x, t), x))^2/w(x, t)^2)+a[10]*(diff(w(x, t), x))/w(x, t)+a[0]

 (5)
> 

K1 := normal(eval(pde, K))

(w(x, t)^4*(diff(diff(w(x, t), x), x))*a[0]*a[10]+w(x, t)^4*(diff(diff(diff(w(x, t), x), x), x))*a[0]*a[20]+w(x, t)^4*(diff(diff(diff(diff(diff(w(x, t), x), x), x), x), x))*delta*a[20]+w(x, t)^4*(diff(diff(diff(diff(w(x, t), x), x), x), x))*delta*a[10]-3*w(x, t)^3*(diff(diff(w(x, t), x), x))^2*delta*a[10]+w(x, t)^3*(diff(diff(w(x, t), x), x))^2*a[10]*a[20]+w(x, t)^3*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))*a[10]^2+w(x, t)^3*(diff(diff(w(x, t), x), x))*(diff(diff(diff(w(x, t), x), x), x))*a[20]^2-w(x, t)^3*(diff(w(x, t), x))^2*a[0]*a[10]-3*w(x, t)^2*(diff(diff(w(x, t), x), x))^2*(diff(w(x, t), x))*a[20]^2+2*w(x, t)^2*(diff(w(x, t), x))^3*a[0]*a[20]-w(x, t)^2*(diff(w(x, t), x))^2*(diff(diff(diff(w(x, t), x), x), x))*a[20]^2+5*w(x, t)*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))^3*a[20]^2-6*w(x, t)*(diff(w(x, t), x))^4*delta*a[10]+3*w(x, t)*(diff(w(x, t), x))^4*a[10]*a[20]-w(x, t)^3*(diff(diff(w(x, t), x), x))*(diff(w(x, t), t))*a[20]-a[10]*(diff(w(x, t), x))*(diff(w(x, t), t))*w(x, t)^3-2*w(x, t)^3*(diff(w(x, t), x))*(diff(diff(w(x, t), t), x))*a[20]+2*w(x, t)^2*(diff(w(x, t), t))*(diff(w(x, t), x))^2*a[20]-3*w(x, t)^3*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))*a[0]*a[20]-10*w(x, t)^3*(diff(diff(w(x, t), x), x))*(diff(diff(diff(w(x, t), x), x), x))*delta*a[20]-4*w(x, t)^3*(diff(w(x, t), x))*(diff(diff(diff(w(x, t), x), x), x))*delta*a[10]+w(x, t)^3*(diff(w(x, t), x))*(diff(diff(diff(w(x, t), x), x), x))*a[10]*a[20]-5*w(x, t)^3*(diff(w(x, t), x))*(diff(diff(diff(diff(w(x, t), x), x), x), x))*delta*a[20]+30*w(x, t)^2*(diff(diff(w(x, t), x), x))^2*(diff(w(x, t), x))*delta*a[20]+12*w(x, t)^2*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))^2*delta*a[10]-5*w(x, t)^2*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))^2*a[10]*a[20]+20*w(x, t)^2*(diff(w(x, t), x))^2*(diff(diff(diff(w(x, t), x), x), x))*delta*a[20]-60*w(x, t)*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))^3*delta*a[20]-w(x, t)^2*(diff(w(x, t), x))^3*a[10]^2+24*(diff(w(x, t), x))^5*delta*a[20]+(diff(diff(diff(w(x, t), t), x), x))*w(x, t)^4*a[20]+a[10]*(diff(diff(w(x, t), t), x))*w(x, t)^4-2*(diff(w(x, t), x))^5*a[20]^2)/w(x, t)^5 = 0

 (6)
> 

K2 := expand(%)

(diff(diff(w(x, t), x), x))*a[0]*a[10]/w(x, t)+(diff(diff(diff(w(x, t), x), x), x))*a[0]*a[20]/w(x, t)+(diff(diff(diff(diff(diff(w(x, t), x), x), x), x), x))*delta*a[20]/w(x, t)+(diff(diff(diff(diff(w(x, t), x), x), x), x))*delta*a[10]/w(x, t)-3*(diff(diff(w(x, t), x), x))^2*delta*a[10]/w(x, t)^2+(diff(diff(w(x, t), x), x))^2*a[10]*a[20]/w(x, t)^2+(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))*a[10]^2/w(x, t)^2+(diff(diff(w(x, t), x), x))*(diff(diff(diff(w(x, t), x), x), x))*a[20]^2/w(x, t)^2-(diff(w(x, t), x))^2*a[0]*a[10]/w(x, t)^2-3*(diff(diff(w(x, t), x), x))^2*(diff(w(x, t), x))*a[20]^2/w(x, t)^3+2*(diff(w(x, t), x))^3*a[0]*a[20]/w(x, t)^3-(diff(w(x, t), x))^2*(diff(diff(diff(w(x, t), x), x), x))*a[20]^2/w(x, t)^3+5*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))^3*a[20]^2/w(x, t)^4-6*(diff(w(x, t), x))^4*delta*a[10]/w(x, t)^4+3*(diff(w(x, t), x))^4*a[10]*a[20]/w(x, t)^4-(diff(diff(w(x, t), x), x))*(diff(w(x, t), t))*a[20]/w(x, t)^2-a[10]*(diff(w(x, t), x))*(diff(w(x, t), t))/w(x, t)^2-2*(diff(w(x, t), x))*(diff(diff(w(x, t), t), x))*a[20]/w(x, t)^2+2*(diff(w(x, t), t))*(diff(w(x, t), x))^2*a[20]/w(x, t)^3+24*(diff(w(x, t), x))^5*delta*a[20]/w(x, t)^5+20*(diff(w(x, t), x))^2*(diff(diff(diff(w(x, t), x), x), x))*delta*a[20]/w(x, t)^3-60*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))^3*delta*a[20]/w(x, t)^4-3*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))*a[0]*a[20]/w(x, t)^2-10*(diff(diff(w(x, t), x), x))*(diff(diff(diff(w(x, t), x), x), x))*delta*a[20]/w(x, t)^2-4*(diff(w(x, t), x))*(diff(diff(diff(w(x, t), x), x), x))*delta*a[10]/w(x, t)^2+(diff(w(x, t), x))*(diff(diff(diff(w(x, t), x), x), x))*a[10]*a[20]/w(x, t)^2-5*(diff(w(x, t), x))*(diff(diff(diff(diff(w(x, t), x), x), x), x))*delta*a[20]/w(x, t)^2+30*(diff(diff(w(x, t), x), x))^2*(diff(w(x, t), x))*delta*a[20]/w(x, t)^3+12*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))^2*delta*a[10]/w(x, t)^3-5*(diff(diff(w(x, t), x), x))*(diff(w(x, t), x))^2*a[10]*a[20]/w(x, t)^3-(diff(w(x, t), x))^3*a[10]^2/w(x, t)^3+(diff(diff(diff(w(x, t), t), x), x))*a[20]/w(x, t)+a[10]*(diff(diff(w(x, t), t), x))/w(x, t)-2*(diff(w(x, t), x))^5*a[20]^2/w(x, t)^5 = 0

 (7)

NULL

Download series-finding.mw

Dealing with Elliptic Function Output...

Asked by: JoyDivisionMan 30
piecewise pdf elliptic-integral
February 15 2025
0  2

I have constructed two Random Variables (X1 and X2).  They work the way I want.  When I take the abs(X1-X2) and ask for the PDF, I get output with Elliptic functions that seems to be in complex space. I read about the Elliptic functions, but got nowhere -- above my head.

I only want the "real" output on the 0 <= t <= 2 domain.  Does anyone know how to parse this output so that I can see a PDF as a function of t?

> 

NULL
> 

with(Statistics)
> 

fx1 := piecewise(t < -1, 0, And(t >= -1, t <= 1), 2*sqrt(-t^2+1)/Pi, t > 1, 0)
> 

fx2 := piecewise(t < -1, 0, And(t >= -1, t <= 1), 2*sqrt(-t^2+1)/Pi, t > 1, 0)
> 

X1 := RandomVariable(Distribution(PDF = unapply(fx1, t)))
> 

X2 := RandomVariable(Distribution(PDF = unapply(fx2, t)))
> 

Dist := abs(X1-X2)
> 

simplify(PDF(Dist, t))

piecewise(t <= 0, 0, t < 2, (1/3)*((-2*t^3-8*t)*EllipticE(I*t/((2+t)^(1/2)*(-t+2)^(1/2)), I*(-t^2+4)^(1/2)/t)+16*EllipticF(I*t/((2+t)^(1/2)*(-t+2)^(1/2)), I*(-t^2+4)^(1/2)/t)*t+(2*t^3+8*t)*EllipticE(I*(-t^2+4)^(1/2)/t)-16*EllipticK(I*(-t^2+4)^(1/2)/t)*t+infinity)/Pi^2, 2 <= t, 0)

 (1)
 

NULL

Download AbsDiff.mw

What is Idea behind finding parameter?...

Asked by: salim-barzani 655
parameters exponential pde
February 15 2025
1  10

I need to find parameter a[12] any one have any vision for finding parameter  , in p2a must contain 3 exponential but we recieve 19 of them which is something i think it is trail function but trail is give me result so must be a way for finding parameter 

a[12]-pde.mw

How to check solutions to a first-order system of ...

Asked by: zenterix 400
vector substitution equality
February 14 2025
0  2

The document below is long because I go through the context of my question step-by-step. 

The context is solving a first-order system of differential equations where we have complex eigenvalues.

But essentially, my question has to do with the end of this reasoning. I find four solutions x1, x2, x3, and x4, and I would like to check whether they are indeed solutions by subbing them into the original system x'=Ax.

When we make such a substitution, we get an equation of the form yl=yr, where yl and yr are 4x1 vectors. If the equality is true then we can say we have a solution.

My question is about how to check for this equality. I amusing Maple's "Equals" but this doesn't seem to work for all four solutions.

with(LinearAlgebra)

NULL

I explain the entire reasoning that leads to the final vectors that I compare

NULL

The following matrix has complex eigenvalues 3+`&+-`(4*I), each with multiplicity 2.

 

A := `<,>`(`<|>`(3, -4, 1, 0), `<|>`(4, 3, 0, 1), `<|>`(0, 0, 3, -4), `<|>`(0, 0, 4, 3))

Matrix(%id = 36893488151974184588)

 (1)

NULL

Consider the eigenvalue `&lambda;__1` = 3-4*I. The associated eigenvectors are

evec := LinearSolve(A-(3-4*I)*IdentityMatrix(4), `<,>`(0, 0, 0, 0), free = c)

Vector[column](%id = 36893488151974164828)

 (2)

We see this gives us only one linearly independent eigenvector.

 

We can find another one by forming a rank 2 generalized eigenvector from one of the eigenvectors above.

 

The pair of such eigenvectors is a chain of generalized eigenvectors of length 2 and they satisfy

 

(A-I*`&lambda;__1`)*w__2 = w__1*(A-I*`&lambda;__1`)*w__1 and w__1*(A-I*`&lambda;__1`)*w__1 = 0

NULL

These equations come from trying the solution w__1*e^(`&lambda;__1`*t) and (t*w__1+w__2)*exp(`&lambda;__1`*t), where w__1 is an eigenvector for `&lambda;__1`.

The equations above lead to

 

(A-I*`&lambda;__1`)^2*w__2 = 0

 

which we can solve for the generalized eigenvector w__2 as follows

NULL

`w__2,gen` := LinearSolve((A-(3-4*I)*IdentityMatrix(4))^2, `<,>`(0, 0, 0, 0), free = c)

Vector[column](%id = 36893488151974138212)

 (3)

One such generalized eigenvector is

 

w__2 := subs({c[2] = I, c[4] = I}, `w__2,gen`)

Vector[column](%id = 36893488151982985812)

 (4)

which gives us

 

w__1 := (A-(3-4*I)*IdentityMatrix(4)).w__2

Vector[column](%id = 36893488151982979180)

 (5)

which is indeed an eigenvector for `&lambda;__1`.

Thus, we have the two complex solutions

 

s__1 := w__1*exp((3-4*I)*t)

Vector[column](%id = 36893488151982971708)

 (6)

s__2 := (t*w__1+w__2)*exp((3-4*I)*t)

Vector[column](%id = 36893488151982963516)

 (7)

But we want real solutions.

NULL

x__1 := `assuming`([Re(s__1)], [t > 0])

Vector[column](%id = 36893488151974145788)

 (8)

x__2 := `assuming`([Im(s__1)], [t > 0])

Vector[column](%id = 36893488151982991948)

 (9)

x__3 := `assuming`([Re(s__2)], [t > 0])

Vector[column](%id = 36893488151982977980)

 (10)

x__4 := `assuming`([Im(s__2)], [t > 0])

Vector[column](%id = 36893488151982967252)

 (11)

Let's check that each of the vectors above (the real and imaginary parts of our complex solutions) are real solutions.

 

I would like to substitute, say, x__1 into the first order system of differential equations diff(x, x) = Ax.

 

At first I tried simply

NULL

subs({x = x__1}, diff(x, t) = A.x)

0 = Matrix(%id = 36893488151974184588).Vector[column](%id = 36893488151982948468)

 (12)

This doesn't work because the diff(x, t) command evaluates to zero before we do the desired substitution.

 

Then I tried

 

subs({x = x__1}, 'diff(x, t)' = A.x)

diff(Vector[column](%id = 36893488151982948468), t) = Matrix(%id = 36893488151974184588).Vector[column](%id = 36893488151982948468)

 (13)

 

   

result := eval(subs({x = x__1}, 'diff(x, t)' = A.x))

Vector[column](%id = 36893488151991655348) = Vector[column](%id = 36893488151982967612)

 (14)

NULL

NULL

By inspection, it seems the two sides are equal and so x__1 indeed satisfies the system.

 

But how do I get maple to tell me this?

 

evalb(result)

false

 (15)

Equal(rhs(result), lhs(result))

true

 (16)

Okay, but does this use of Equals always work?

 

sols := [x__1, x__2, x__3, x__4]; for i to 4 do result := eval(subs({x = sols[i]}, 'diff(x, t)' = A.x)); print(simplify(result))*print(Equal(rhs(result), lhs(result))) end do

true

 (17)

 

Apparently not.``

 

What is going on here?

It sure looks like the righthand and lefthand sides of each expression are the same.

Download evs_equal.mw

How to create sections of a document that can be e...

Asked by: zenterix 400
document editing section
February 14 2025
0  2

I am creating a Maple Document, and I would like to separate the document into sections.

One way that looks like what I could use is what I see in parts of the documentation, as below

So here the document looks like an outline.

How is this actually done?

Are there other ways of doing something like this?

solve issues

by: Yehor Volovodenko 10 Maple
solve parameters
February 14 2025
2  4

Greetings, dear Maple developers!

I, Yegor Volovodenko, together with my supervisor, Igor Zinoviev, Associate Professor, Head of the Department of General Mathematics at ZNU, am conducting research on ‘Solving equations and inequalities of elementary mathematics by means of computer mathematics: opportunities, problems and ways to solve them’. I express my sincere gratitude for your work on the development and popularisation of mathematics, for the constant improvement and enhancement of CAS in particular.

Unfortunately, during the study, we found some, in our opinion, shortcomings in the work of Maxima algorithms for solving elementary equations and inequalities

  1. Solving equations with parameters;



    In these two cases, the solution obtained is formal, without analysing the cases when the coefficient of a variable is zero. Thus, not all solutions of equations with a parameter are obtained. I believe that if the user is not familiar with the methodology for solving equations with parameters, he or she may lose some solutions that may be important for further work. So the solution can be considered incomplete.

 

  1. Solving equations with two variables.
     

In the equation with two variables, the answer was given in complex numbers, which is incomprehensible to a person who is not familiar with complex numbers. There are also comments on the course of the solution, in this equation it was necessary to select the square of the difference, and then solve the equation x^2+(y-2)^2=0, getting x=0, y=2.

I hope that the results I have obtained (the identified shortcomings) will help to correct the work of the algorithms and improve the work of the Maple system.

Sincerely, Yegor Volovodenko, Igor Zinoviev.

How to use fsolve with numerical solutions of ODEs...

Asked by: fnavarro 45
ode fsolve numeric
February 14 2025
1  2

Hi! I have a problem that is driving me crazy. It is related to shooting methods for ODEs with Maple. Essentially what one needs is to find the solution to two equations simultaneously, defined in terms of the numerical solution to a second-order differential equation. I have used fsolve and also DirectSearch but nothing works. Graphically one can see the solution to the problem exists.

I include the Maple file. Thank you. 

Example of shooting

> 

restart
> 

with(plots)

ODE to be solved

> 

eq_f_bar_compact := (1-x)^4*((diff(f_bar(x), x, x))*x^2+4*(diff(f_bar(x), x))*x+2*f_bar(x))-2*(1-x)^3*((diff(f_bar(x), x))*x^2+2*f_bar(x)*x)-(-(4*((diff(f_bar(x), x))*x^2+2*f_bar(x)*x))*(1-x)^5+(1/5)*sin(f_bar(x)*x^2)*x^3+4*cos(f_bar(x)*x^2)^2*((diff(f_bar(x), x))*x^2+2*f_bar(x)*x)*(1-x)^5-4*cos(f_bar(x)*x^2)*sin(f_bar(x)*x^2)*(1-x)^2*x+4*cos(f_bar(x)*x^2)*sin(f_bar(x)*x^2)*((diff(f_bar(x), x))*x^2+2*f_bar(x)*x)^2*(1-x)^6*x+((diff(f_bar(x), x))*x^2+2*f_bar(x)*x)*(1-x)^3*x^2+(1/5)*cos(f_bar(x)*x^2)*sin(f_bar(x)*x^2)*x^3)/(x*(4*cos(f_bar(x)*x^2)^2*(1-x)^2-4*(1-x)^2-x^2))

(1-x)^4*((diff(diff(f_bar(x), x), x))*x^2+4*(diff(f_bar(x), x))*x+2*f_bar(x))-2*(1-x)^3*((diff(f_bar(x), x))*x^2+2*f_bar(x)*x)-(-4*((diff(f_bar(x), x))*x^2+2*f_bar(x)*x)*(1-x)^5+(1/5)*sin(f_bar(x)*x^2)*x^3+4*cos(f_bar(x)*x^2)^2*((diff(f_bar(x), x))*x^2+2*f_bar(x)*x)*(1-x)^5-4*cos(f_bar(x)*x^2)*sin(f_bar(x)*x^2)*(1-x)^2*x+4*cos(f_bar(x)*x^2)*sin(f_bar(x)*x^2)*((diff(f_bar(x), x))*x^2+2*f_bar(x)*x)^2*(1-x)^6*x+((diff(f_bar(x), x))*x^2+2*f_bar(x)*x)*(1-x)^3*x^2+(1/5)*cos(f_bar(x)*x^2)*sin(f_bar(x)*x^2)*x^3)/(x*(4*cos(f_bar(x)*x^2)^2*(1-x)^2-4*(1-x)^2-x^2))

 (1)
> 

indets(eq_f_bar_compact)

{x, cos(f_bar(x)*x^2), diff(diff(f_bar(x), x), x), diff(f_bar(x), x), f_bar(x), sin(f_bar(x)*x^2)}

 (2)

Range x in [0,1]

Matching point

> 

x_match := .5

.5

 (3)

Let us impose the BCs at a point a little larger than 0 for x=0 and a little smaller than 1 for x=1

> 

eps := 10^(-3)

1/1000

 (4)

*********************
Shooting from x=0
*********************

Expansion at x=0

> 

origin_expansion_for_f_bar := {f_bar(x) = f0+2*f0*x+(-(4/3)*f0^3+(89/30)*f0)*x^2+(-(16/3)*f0^3+(58/15)*f0)*x^3+((14/3)*f0^5-(797/60)*f0^3+(11201/2400)*f0)*x^4+(28*f0^5-(791/30)*f0^3+(6403/1200)*f0)*x^5+(-(328/15)*f0^7+(1468/15)*f0^5-(410533/9000)*f0^3+(2103149/360000)*f0)*x^6}

{f_bar(x) = f0+2*f0*x+(-(4/3)*f0^3+(89/30)*f0)*x^2+(-(16/3)*f0^3+(58/15)*f0)*x^3+((14/3)*f0^5-(797/60)*f0^3+(11201/2400)*f0)*x^4+(28*f0^5-(791/30)*f0^3+(6403/1200)*f0)*x^5+(-(328/15)*f0^7+(1468/15)*f0^5-(410533/9000)*f0^3+(2103149/360000)*f0)*x^6}

 (5)

Boundary condition at x=eps. c0 is a free parameter (to be determined)

> 

BCs_x_0 := subs(x = eps, {f_bar(x) = subs({f0 = c0}, subs(origin_expansion_for_f_bar, f_bar(x))), (D(f_bar))(x) = subs({f0 = c0}, eval(diff(subs(origin_expansion_for_f_bar, f_bar(x)), x)))})
> 

dsys_x_0 := `union`({eq_f_bar_compact}, BCs_x_0)

Numerical solution of the equation

> 

fun_sol_x_0 := dsolve(dsys_x_0, numeric, parameters = [c0])

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 := [c0 = c0]; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 28, [( 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) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 1, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (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.10e-2, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = 0.10e-2, (6) = .0, (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..3, {(1) = 1.00200297053801*c0-0.133867997641228e-5*c0^3+0.469476453333333e-11*c0^5-0.218666666666667e-16*c0^7, (2) = 2.00594495202838*c0-0.268271993210702e-2*c0^3+0.188072538666667e-7*c0^5-0.131200000000000e-12*c0^7, (3) = Float(undefined)})), ( 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..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8]), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..2, {(1, 1) = .0, (1, 2) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = f_bar(x), Y[2] = diff(f_bar(x),x)]`; YP[2] := -((1-X)^4*(4*X*Y[2]+2*Y[1])-2*(1-X)^3*(X^2*Y[2]+2*X*Y[1])-(-4*(X^2*Y[2]+2*X*Y[1])*(1-X)^5+(1/5)*sin(Y[1]*X^2)*X^3+4*cos(Y[1]*X^2)^2*(X^2*Y[2]+2*X*Y[1])*(1-X)^5-4*cos(Y[1]*X^2)*sin(Y[1]*X^2)*(1-X)^2*X+4*cos(Y[1]*X^2)*sin(Y[1]*X^2)*(X^2*Y[2]+2*X*Y[1])^2*(1-X)^6*X+(X^2*Y[2]+2*X*Y[1])*(1-X)^3*X^2+(1/5)*cos(Y[1]*X^2)*sin(Y[1]*X^2)*X^3)/(X*(4*cos(Y[1]*X^2)^2*(1-X)^2-4*(1-X)^2-X^2)))/((1-X)^4*X^2); YP[1] := Y[2]; 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] = f_bar(x), Y[2] = diff(f_bar(x),x)]`; YP[2] := -((1-X)^4*(4*X*Y[2]+2*Y[1])-2*(1-X)^3*(X^2*Y[2]+2*X*Y[1])-(-4*(X^2*Y[2]+2*X*Y[1])*(1-X)^5+(1/5)*sin(Y[1]*X^2)*X^3+4*cos(Y[1]*X^2)^2*(X^2*Y[2]+2*X*Y[1])*(1-X)^5-4*cos(Y[1]*X^2)*sin(Y[1]*X^2)*(1-X)^2*X+4*cos(Y[1]*X^2)*sin(Y[1]*X^2)*(X^2*Y[2]+2*X*Y[1])^2*(1-X)^6*X+(X^2*Y[2]+2*X*Y[1])*(1-X)^3*X^2+(1/5)*cos(Y[1]*X^2)*sin(Y[1]*X^2)*X^3)/(X*(4*cos(Y[1]*X^2)^2*(1-X)^2-4*(1-X)^2-X^2)))/((1-X)^4*X^2); YP[1] := Y[2]; 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), ( 28 ) = (0) ] )) ] ); _y0 := Array(0..3, {(1) = 0.100000000000000e-2, (2) = 1.00200297053801*c0-0.133867997641228e-5*c0^3+0.469476453333333e-11*c0^5-0.218666666666667e-16*c0^7, (3) = 2.00594495202838*c0-0.268271993210702e-2*c0^3+0.188072538666667e-7*c0^5-0.131200000000000e-12*c0^7}); _vmap := array( 1 .. 2, [( 1 ) = (1), ( 2 ) = (2) ] ); _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 elif type(_xin, `=`) and lhs(_xin) = "setdatacallback" then if not type(rhs(_xin), 'nonegint') then error "data callback must be a nonnegative integer (address)" end if; _dtbl[1][28] := rhs(_xin) 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 <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> 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 <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> 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) = [x, f_bar(x), diff(f_bar(x), x)], (4) = [c0 = c0]}); _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

 (6)
> 

sol_x_0 := proc (c) fun_sol_x_0(parameters = [c]); fun_sol_x_0(x_match) end proc

proc (c) fun_sol_x_0(parameters = [c]); fun_sol_x_0(x_match) end proc

 (7)
> 

sol_x_0(-.5)

[x = .5, f_bar(x) = HFloat(-1.5657446902999976), diff(f_bar(x), x) = HFloat(-4.035127384355003)]

 (8)
> 

plot('subs(sol_x_0(c), f_bar(x))', c = -.7 .. -.1)

 

*********************
Shooting from x=1
*********************

> 

asymp_expansion_for_f_bar := {f_bar(x) = -Pi+2*Pi*(x-1)+f2*(x-1)^2+(-8*Pi-4*f2)*(x-1)^3+((1/120)*f2^3+3*f2^2*Pi*(1/40)+(1/120)*(27*Pi^2+1200)*f2+9*Pi^3*(1/40)+25*Pi)*(x-1)^4+(-(1/20)*f2^3-9*f2^2*Pi*(1/20)+(1/20)*(-27*Pi^2-400)*f2-27*Pi^3*(1/20)-54*Pi)*(x-1)^5+(243*Pi^5*(1/12800)+81*Pi^4*f2*(1/2560)+(1/1280)*(9*(3*f2^2+352))*Pi^3+(1/1280)*(9*(f2^3+352*f2))*Pi^2+(1/2560)*(3*f2^4+2112*f2^2+250880)*Pi+(1/12800)*f2^5+11*f2^3*(1/120)+35*f2)*(x-1)^6}

{f_bar(x) = -Pi+2*Pi*(x-1)+f2*(x-1)^2+(-8*Pi-4*f2)*(x-1)^3+((1/120)*f2^3+(3/40)*f2^2*Pi+(1/120)*(27*Pi^2+1200)*f2+(9/40)*Pi^3+25*Pi)*(x-1)^4+(-(1/20)*f2^3-(9/20)*f2^2*Pi+(1/20)*(-27*Pi^2-400)*f2-(27/20)*Pi^3-54*Pi)*(x-1)^5+((243/12800)*Pi^5+(81/2560)*Pi^4*f2+(9/1280)*(3*f2^2+352)*Pi^3+(9/1280)*(f2^3+352*f2)*Pi^2+(1/2560)*(3*f2^4+2112*f2^2+250880)*Pi+(1/12800)*f2^5+(11/120)*f2^3+35*f2)*(x-1)^6}

 (9)

Boundary condition at x=1-eps. c1 is a free parameter (to be determined)

> 

BCs_x_1 := subs(x = 1-eps, {f_bar(x) = subs({f2 = c1}, subs(asymp_expansion_for_f_bar, f_bar(x))), (D(f_bar))(x) = subs({f2 = c1}, eval(diff(subs(asymp_expansion_for_f_bar, f_bar(x)), x)))})
> 

dsys_x_1 := `union`({eq_f_bar_compact}, BCs_x_1)

Numerical solution of the equation

> 

fun_sol_x_1 := dsolve(dsys_x_1, numeric, parameters = [c1])

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 := [c1 = c1]; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 28, [( 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) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 1, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (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) = .999, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .999, (6) = .0, (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..3, {(1) = -3.14787581367851+0.100401225404747e-5*c1+0.838349439565595e-14*c1^3+0.237036411565938e-12*c1^2+0.368155389092554e-20*c1^4+0.781250000000000e-22*c1^5, (2) = 6.28310956583112-0.201204904963884e-2*c1-0.335842997072690e-10*c1^3-0.949565854663041e-9*c1^2-0.220893233455532e-16*c1^4-0.468750000000000e-18*c1^5, (3) = Float(undefined)})), ( 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..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8]), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8])]), ( 8 ) = ([Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..2, {(1, 1) = .0, (1, 2) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = f_bar(x), Y[2] = diff(f_bar(x),x)]`; YP[2] := -((1-X)^4*(4*X*Y[2]+2*Y[1])-2*(1-X)^3*(X^2*Y[2]+2*X*Y[1])-(-4*(X^2*Y[2]+2*X*Y[1])*(1-X)^5+(1/5)*sin(Y[1]*X^2)*X^3+4*cos(Y[1]*X^2)^2*(X^2*Y[2]+2*X*Y[1])*(1-X)^5-4*cos(Y[1]*X^2)*sin(Y[1]*X^2)*(1-X)^2*X+4*cos(Y[1]*X^2)*sin(Y[1]*X^2)*(X^2*Y[2]+2*X*Y[1])^2*(1-X)^6*X+(X^2*Y[2]+2*X*Y[1])*(1-X)^3*X^2+(1/5)*cos(Y[1]*X^2)*sin(Y[1]*X^2)*X^3)/(X*(4*cos(Y[1]*X^2)^2*(1-X)^2-4*(1-X)^2-X^2)))/((1-X)^4*X^2); YP[1] := Y[2]; 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] = f_bar(x), Y[2] = diff(f_bar(x),x)]`; YP[2] := -((1-X)^4*(4*X*Y[2]+2*Y[1])-2*(1-X)^3*(X^2*Y[2]+2*X*Y[1])-(-4*(X^2*Y[2]+2*X*Y[1])*(1-X)^5+(1/5)*sin(Y[1]*X^2)*X^3+4*cos(Y[1]*X^2)^2*(X^2*Y[2]+2*X*Y[1])*(1-X)^5-4*cos(Y[1]*X^2)*sin(Y[1]*X^2)*(1-X)^2*X+4*cos(Y[1]*X^2)*sin(Y[1]*X^2)*(X^2*Y[2]+2*X*Y[1])^2*(1-X)^6*X+(X^2*Y[2]+2*X*Y[1])*(1-X)^3*X^2+(1/5)*cos(Y[1]*X^2)*sin(Y[1]*X^2)*X^3)/(X*(4*cos(Y[1]*X^2)^2*(1-X)^2-4*(1-X)^2-X^2)))/((1-X)^4*X^2); YP[1] := Y[2]; 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), ( 28 ) = (0) ] )) ] ); _y0 := Array(0..3, {(1) = .999000000000000, (2) = -3.14787581367851+0.100401225404747e-5*c1+0.838349439565595e-14*c1^3+0.237036411565938e-12*c1^2+0.368155389092554e-20*c1^4+0.781250000000000e-22*c1^5, (3) = 6.28310956583112-0.201204904963884e-2*c1-0.335842997072690e-10*c1^3-0.949565854663041e-9*c1^2-0.220893233455532e-16*c1^4-0.468750000000000e-18*c1^5}); _vmap := array( 1 .. 2, [( 1 ) = (1), ( 2 ) = (2) ] ); _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 elif type(_xin, `=`) and lhs(_xin) = "setdatacallback" then if not type(rhs(_xin), 'nonegint') then error "data callback must be a nonnegative integer (address)" end if; _dtbl[1][28] := rhs(_xin) 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 <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> 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 <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> 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) = [x, f_bar(x), diff(f_bar(x), x)], (4) = [c1 = c1]}); _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

 (10)
> 

sol_x_1 := proc (c) fun_sol_x_1(parameters = [c]); fun_sol_x_1(x_match) end proc

proc (c) fun_sol_x_1(parameters = [c]); fun_sol_x_1(x_match) end proc

 (11)
> 

sol_x_1(.8)

[x = .5, f_bar(x) = HFloat(-2.9927634069258033), diff(f_bar(x), x) = HFloat(2.802868299574931)]

 (12)
> 

plot('subs(sol_x_1(c), f_bar(x))', c = .1 .. .8)

 

Now we generate twi numerical functions in two variables (c0,c1) to get the solution: we have to find the zeroes of both functions simultaneously

> 

zero_function1 := proc (c0_, c1_) subs(sol_x_0(c0_), f_bar(x))-subs(sol_x_1(c1_), f_bar(x)) end proc
> 

zero_function1(-.5, .1)

HFloat(1.7905971356833628)

 (13)
> 

NULL
> 

zero_function2 := proc (c0_, c1_) subs(sol_x_0(c0_), diff(f_bar(x), x))-subs(sol_x_1(c1_), diff(f_bar(x), x)) end proc
> 

NULL
> 

zero_function2(-.5, .1)

HFloat(-8.481253652259962)

 (14)

We first try fsolve

> 

fsolve({'zero_function1(a, b)', 'zero_function2(a, b)'}, {a = -.6, b = 3})

fsolve({zero_function1(a, b), zero_function2(a, b)}, {a = -.6, b = 3})

 (15)

It does not work

But it is clear a solution EXISTS!

> 

plot1 := plot3d('zero_function1(a, b)', a = -.62 .. -.6, b = 2.85 .. 2.95, color = red)

 
> 

plot2 := plot3d('zero_function2(a, b)', a = -.62 .. -.6, b = 2.85 .. 2.95, color = blue)

 
> 

plot3 := plot3d(0, a = -.62 .. -.6, b = 2.85 .. 2.95, color = green)

 
> 

display({plot1, plot2, plot3})

 

There is a solution!

Let us try the package DirectSearch

> 

with(DirectSearch)

[BoundedObjective, CompromiseProgramming, DataFit, ExponentialWeightedSum, GlobalOptima, GlobalSearch, Minimax, ModifiedTchebycheff, Search, SolveEquations, WeightedProduct, WeightedSum]

 (16)

Let us try to find it as a minimum

> 

GlobalSearch('zero_function1(a, b)^2+zero_function2(a, b)^2', pointrange = [a = -.7 .. .5, b = 2 .. 4])

Error, (in dsolve/numeric/process_parameters) 'parameters' must be specified as a list of numeric values

  

We get a mistake.

Any idea how to solve the error (and the problem of finding the solution to "{zero_function1(a,b)=0, ""zero_function2(a,b)=0}"?

Download shooting_question.mw

How can we collect coefficients of different power...

Asked by: Hamza Rafiq 30
homework collect differential-equations
February 14 2025
0  2

How can we collect coefficients of different powers of involved variables form Eq, to get system of algebraic equations, like x,y,z,t but these are also appear in JacobiCN function? Maple file link is attached below, please help me to solve this issue. 

collection_of_coefficients.mw

How can check pdetest when function is not short? ...

Asked by: salim-barzani 655
explore pde differential-equations
February 14 2025
1  4

the function is true but i want to be sure when i use pdetest must give me zero, but there must be a way for checking such function, please if your pc not strong don't click the command pdetest, i want use explore for such function but i am not sure it work or not, becuase the graph are a little bit strange  and long , i want  a way for easy plotting and visualization of such graph , can anyone help for solve this issue?

 sol.mw

How to solve a first order ODE numerically?...

Asked by: mary120 75
ode dsolve numeric
February 13 2025
1  16

Hi,

I want to solve a first order ODE and (i) plot phi vs. x and then (ii) export data of the plotted curve as an T-shape ASCII file. I try to do this by Maple, but there are some errors and I couldn’t get the mentioned curve and ASCII file. Please, help me to remove errors:
(dphi/dx)**2+2*V(phi)=0,

where V(phi)= (1-alpha)*M^2 - (1-alpha)*M*sqrt(M^2 - 2*phi) + mu*(mu + beta*nu)*(1 - exp(phi/(mu + beta*nu))) + 
            (nu/beta)*(mu + beta*nu)*(1 - exp(beta*phi/(mu + beta*nu))) + (alpha/gamma1)*(1 - exp(-gamma1*phi)),

where we assume mu := 0.01, nu := 1 - mu, beta := 0.05, alpha := 0.3, M := sqrt(0.704), gamma1 := 0.001.
As it is seen from the attached figure, V(phi) and dV(phi)/dx =0 for phi=0 and phi=phi_m (two extreme points of V(phi)) and d^V(phi)/dphi^2<0 at phi=0 and phi=phi_m (phi and x are real). 

From the attached figure, it seems that phi=phi_m at x=0 and phi->0 as x-> infinity. (Hint: I expect the plot of V(phi) vs. phi and phi vs. x be similar to the curve "A" in the attached file).

Thanks.
rk4.mw

How apply conjugate for getting result of pdetest ...

Asked by: salim-barzani 655
pde differential-equations conjugate
February 13 2025
2  2

i did pdetest without conjugate like the paper did i get zero but when i did pde test  with conjugate i didn't where is my problem 
i will do without conjugate but how change p[2]=conjugate(p[1])

> 

restart
> 

with(PDEtools)
> 

with(LinearAlgebra)
> 

NULL
> 

with(SolveTools)
> 

undeclare(prime)

`There is no more prime differentiation variable; all derivatives will be displayed as indexed functions`

 (1)
> 

declare(u(x, y, z, t))

u(x, y, z, t)*`will now be displayed as`*u

 (2)
> 

declare(f(x, y, z, t))

f(x, y, z, t)*`will now be displayed as`*f

 (3)
> 

pde := -4*(diff(u(x, y, z, t), x, t))+diff(u(x, y, z, t), `$`(x, 3), z)+3*alpha*(diff(u(x, y, z, t), `$`(y, 2)))+4*(diff(u(x, y, z, t), x))*(diff(u(x, y, z, t), x, z))+2*(diff(u(x, y, z, t), `$`(x, 2)))*(diff(u(x, y, z, t), z))

-4*(diff(diff(u(x, y, z, t), t), x))+diff(diff(diff(diff(u(x, y, z, t), x), x), x), z)+3*alpha*(diff(diff(u(x, y, z, t), y), y))+4*(diff(u(x, y, z, t), x))*(diff(diff(u(x, y, z, t), x), z))+2*(diff(diff(u(x, y, z, t), x), x))*(diff(u(x, y, z, t), z))

 (4)
> 

pde_nonlinear, pde_linear := selectremove(proc (term) options operator, arrow; has((eval(term, u(x, y, z, t) = a*u(x, y, z, t)))/a, a) end proc, pde)

4*(diff(u(x, y, z, t), x))*(diff(diff(u(x, y, z, t), x), z))+2*(diff(diff(u(x, y, z, t), x), x))*(diff(u(x, y, z, t), z)), -4*(diff(diff(u(x, y, z, t), t), x))+diff(diff(diff(diff(u(x, y, z, t), x), x), x), z)+3*alpha*(diff(diff(u(x, y, z, t), y), y))

 (5)
> 

thetai := t*w[i]+y*p[i]+x+z

t*w[i]+y*p[i]+x+z

 (6)
> 

eqw := w[i] = 3*alpha*p[i]^2*(1/4)

w[i] = (3/4)*alpha*p[i]^2

 (7)
> 

Bij := proc (i, j) options operator, arrow; 4/((p[i]-p[j])^2*alpha) end proc

proc (i, j) options operator, arrow; 4/((p[i]-p[j])^2*alpha) end proc

 (8)
> 

theta1 := normal(eval(eval(thetai, eqw), i = 1)); theta2 := normal(eval(eval(thetai, eqw), i = 2))

(3/4)*alpha*t*p[1]^2+y*p[1]+x+z

  

(3/4)*alpha*t*p[2]^2+y*p[2]+x+z

 (9)
> 

eqf := f(x, y, z, t) = theta1*theta2+4/((p[1]-p[2])^2*alpha)

f(x, y, z, t) = ((3/4)*alpha*t*p[1]^2+y*p[1]+x+z)*((3/4)*alpha*t*p[2]^2+y*p[2]+x+z)+4/((p[1]-p[2])^2*alpha)

 (10)
> 

eq17 := u(x, y, z, t) = 2*(diff(ln(f(x, y, z, t)), x))

u(x, y, z, t) = 2*(diff(f(x, y, z, t), x))/f(x, y, z, t)

 (11)
> 

eqt := eval(eq17, eqf)

u(x, y, z, t) = 2*((3/4)*alpha*t*p[2]^2+y*p[2]+2*x+2*z+(3/4)*alpha*t*p[1]^2+y*p[1])/(((3/4)*alpha*t*p[1]^2+y*p[1]+x+z)*((3/4)*alpha*t*p[2]^2+y*p[2]+x+z)+4/((p[1]-p[2])^2*alpha))

 (12)
> 

NULL
> 

pdetest(eqt, pde)

0

 (13)
> 

NULL

Download p1.mw

Plotting an AM waveform...

Asked by: 3cells 10
syntax plotting incomplete-question
February 13 2025
2  4

I am in the process of trying to learn Maple. As an exercise I am trying to plot an amplitude modulated waveform. No problem plotting the waveform of the carrier and the modulating signal but for some reason (user error on my part) I can't plot the combined waveform y given by. 

y = [1 + M*cos(2*Pi*f__m*t)]*A*sin(2*Pi*f__c*t)

Where the following values have been set M = 0.5, A = 1, f_m = 2, f_c = 10.

Plotting y using the following

plot(y, t = 0 .. 1, title = "Graph of Carrier", labels = ["time (seconds)", "Amplitude (Volts)"], color = red);

Gives the following error. Having tried various combinations of variables, values, single parameter functions etc I need a pointer as to where I am going wrong. 

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct
What I am expecting is a graph similar to the below derived in Excel? My version of Maple is 2023.2

Error, (in dsolve/numeric/bvp) initial Newton iter...

Asked by: Saha 80
ode dsolve numeric
February 13 2025
1  3

VC1.mw

How to solve this error?

How do convert language to Matlab when expression...

Asked by: salim-barzani 655
codegeneration matlab
February 13 2025
0  0

i don't know how apply conversation language to matlab in righ hand side  don't show up to do conversation language for short is come up but for this not 

> 

restart
> 

K := (2*(k[1]*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1])+((p[1]-p[2])^2*alpha-(k[1]-k[2])^2)*(k[1]+k[2])*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2])/((p[1]-p[2])^2*alpha-(k[1]+k[2])^2)+k[2]*exp((3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2])+((p[2]-p[3])^2*alpha-(k[2]-k[3])^2)*(k[2]+k[3])*exp((3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/((p[2]-p[3])^2*alpha-(k[2]+k[3])^2)+((p[1]-p[2])^2*alpha-(k[1]-k[2])^2)*((p[1]-p[3])^2*alpha-(k[1]-k[3])^2)*((p[2]-p[3])^2*alpha-(k[2]-k[3])^2)*(k[1]+k[2]+k[3])*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/(((p[1]-p[2])^2*alpha-(k[1]+k[2])^2)*((p[1]-p[3])^2*alpha-(k[1]+k[3])^2)*((p[2]-p[3])^2*alpha-(k[2]+k[3])^2))+((p[1]-p[3])^2*alpha-(k[1]-k[3])^2)*(k[1]+k[3])*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/((p[1]-p[3])^2*alpha-(k[1]+k[3])^2)+k[3]*exp((3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])))/(1+exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1])+((p[1]-p[2])^2*alpha-(k[1]-k[2])^2)*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2])/((p[1]-p[2])^2*alpha-(k[1]+k[2])^2)+exp((3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2])+((p[2]-p[3])^2*alpha-(k[2]-k[3])^2)*exp((3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/((p[2]-p[3])^2*alpha-(k[2]+k[3])^2)+((p[1]-p[2])^2*alpha-(k[1]-k[2])^2)*((p[1]-p[3])^2*alpha-(k[1]-k[3])^2)*((p[2]-p[3])^2*alpha-(k[2]-k[3])^2)*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/(((p[1]-p[2])^2*alpha-(k[1]+k[2])^2)*((p[1]-p[3])^2*alpha-(k[1]+k[3])^2)*((p[2]-p[3])^2*alpha-(k[2]+k[3])^2))+((p[1]-p[3])^2*alpha-(k[1]-k[3])^2)*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/((p[1]-p[3])^2*alpha-(k[1]+k[3])^2)+exp((3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3]))

2*(k[1]*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1])+((p[1]-p[2])^2*alpha-(k[1]-k[2])^2)*(k[1]+k[2])*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2])/((p[1]-p[2])^2*alpha-(k[1]+k[2])^2)+k[2]*exp((3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2])+((p[2]-p[3])^2*alpha-(k[2]-k[3])^2)*(k[2]+k[3])*exp((3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/((p[2]-p[3])^2*alpha-(k[2]+k[3])^2)+((p[1]-p[2])^2*alpha-(k[1]-k[2])^2)*((p[1]-p[3])^2*alpha-(k[1]-k[3])^2)*((p[2]-p[3])^2*alpha-(k[2]-k[3])^2)*(k[1]+k[2]+k[3])*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/(((p[1]-p[2])^2*alpha-(k[1]+k[2])^2)*((p[1]-p[3])^2*alpha-(k[1]+k[3])^2)*((p[2]-p[3])^2*alpha-(k[2]+k[3])^2))+((p[1]-p[3])^2*alpha-(k[1]-k[3])^2)*(k[1]+k[3])*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/((p[1]-p[3])^2*alpha-(k[1]+k[3])^2)+k[3]*exp((3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3]))/(1+exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1])+((p[1]-p[2])^2*alpha-(k[1]-k[2])^2)*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2])/((p[1]-p[2])^2*alpha-(k[1]+k[2])^2)+exp((3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2])+((p[2]-p[3])^2*alpha-(k[2]-k[3])^2)*exp((3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/((p[2]-p[3])^2*alpha-(k[2]+k[3])^2)+((p[1]-p[2])^2*alpha-(k[1]-k[2])^2)*((p[1]-p[3])^2*alpha-(k[1]-k[3])^2)*((p[2]-p[3])^2*alpha-(k[2]-k[3])^2)*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[2]*t*alpha*p[2]^2+(1/4)*t*k[2]^3+k[2]*p[2]*y+k[2]*x+k[2]*z+eta[2]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/(((p[1]-p[2])^2*alpha-(k[1]+k[2])^2)*((p[1]-p[3])^2*alpha-(k[1]+k[3])^2)*((p[2]-p[3])^2*alpha-(k[2]+k[3])^2))+((p[1]-p[3])^2*alpha-(k[1]-k[3])^2)*exp((3/4)*k[1]*t*alpha*p[1]^2+(1/4)*t*k[1]^3+k[1]*p[1]*y+k[1]*x+k[1]*z+eta[1]+(3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3])/((p[1]-p[3])^2*alpha-(k[1]+k[3])^2)+exp((3/4)*k[3]*t*alpha*p[3]^2+(1/4)*t*k[3]^3+k[3]*p[3]*y+k[3]*x+k[3]*z+eta[3]))

 (1)
 

NULL

Download convert-to-matlab.mw

First 7 8 9 10 11 12 13 Last Page 9 of 2191