mmcdara

5974 Reputation

17 Badges

7 years, 265 days

MaplePrimes Activity


These are answers submitted by mmcdara

this proposal.mw ?

restart

 # Parameters always positive

assume(0 < sigma__d, 0 < sigma__v):
interface(showassumed=0):

Diff('lambda__1', sigma__v) = sqrt(5)/(5*sigma__d):
f__1 := rhs(%):
Diff('lambda__1', sigma__d) = -sqrt(5)*sigma__v/(5*sigma__d^2):
f__2 := abs(rhs(%)):

use plots in
  display(
    inequal(
      f__1-f__2 > 0
      , sigma__v=0.001..10, sigma__d=0.001..10
      , color="Chartreuse"
      , optionsexcluded=[color="Niagara DarkOrchid"]
    )
    ,
    # For a more complex case than yours
    # (here a simple plot(solve(f__1=f__2, sigma__d), sigma__d=0.001..10)
    # would be enough)
    implicitplot(
      f__1=f__2
      , sigma__v=0.001..10, sigma__d=0.001..10
      , color=blue
      , thickness=5
      , grid=[100, 100]
      , legend=typeset('f__1' = 'f__2')
    )
    , textplot([3, 7, typeset('f__1' > 'f__2')], font=[Times, bold, 14])
    , textplot([7, 3, typeset('f__1' < 'f__2')], font=[Times, bold, 14], color=white)
    , labels=[typeset(sigma__v), typeset(sigma__d)]
  )
end use;

 

Download proposal.mw

(customize it as you want) 

See also add_on.mw

@jalal 


To get this histogram you present:

  • Given the nature of your data I suppose they are already consist in a serie of elements of the type [C[k], N[k]] where
    • C[k] is the range of the kth class,
    • N[k] is the number of elements C[k] contains.

This is what I call synthetic data by opposition to raw data (see below).
 

  • If it is so you have two possibilities:
    1. To build an Statistics:-Histogram you must have raw data, not synthetic data.
      So either you have these raw data or you don't.
      • If you have them, let R the list/vector of incomes, just use Statistics:-Histogram with the ad hoc customization (first plot in the attached file).
         
      • If you do not have them, you can recreate fake raw data from synthetic data so that their histogram is the one you want (second plot in the attached file).
         
    2. If you only have synthetic data and bo not necessarily want to use Statistics:-Histogram, then you can proceed as explaine in the third plot in the attached file.



 

restart


CASE 1: YOU HAVE RAW DATA


Simulation of hypothetic raw data

R := Statistics:-Sample(Triangular(6, 97, 53), 10^3):


Raw data histogram

C := [seq(0..100, 10)]:
Statistics:-Histogram(
  R
  , frequencyscale=absolute
  , binbounds=C
  , axis[1]=[tickmarks=C]
  , view=[-1..101, default]
  , color="Pink"
  , transparency=0.5
  , labels=["Salaires (milliers de \$)", "Nombre d'employés"]
  , labeldirections=[default, vertical]
  , labelfont=[Tahoma, 10]
  , size=[700, 400]
)

 


CASE 2: YOU DO NOT HAVE RAW DATA BUT SYNTHETIC DATA ONLY

... AND YOU WANT TO USE  Statistics:-Histogram


Here are the synthetic data which correspond to the raw data R above

Classes       := [seq(C[k]..C[k+1], k=1..numelems(C)-1)]:
SyntheticData := Statistics:-TallyInto(R, Classes):
SyntheticData := Matrix(numelems(Classes), 2, (i, j) -> `if`(j=1, lhs(SyntheticData[i]), rhs(SyntheticData[i])));
 

SyntheticData := Matrix(10, 2, {(1, 1) = HFloat(0.0) .. HFloat(10.0), (1, 2) = 3, (2, 1) = HFloat(10.0) .. HFloat(20.0), (2, 2) = 43, (3, 1) = HFloat(20.0) .. HFloat(30.0), (3, 2) = 92, (4, 1) = HFloat(30.0) .. HFloat(40.0), (4, 2) = 141, (5, 1) = HFloat(40.0) .. HFloat(50.0), (5, 2) = 195, (6, 1) = HFloat(50.0) .. HFloat(60.0), (6, 2) = 218, (7, 1) = HFloat(60.0) .. HFloat(70.0), (7, 2) = 152, (8, 1) = HFloat(70.0) .. HFloat(80.0), (8, 2) = 95, (9, 1) = HFloat(80.0) .. HFloat(90.0), (9, 2) = 51, (10, 1) = HFloat(90.0) .. HFloat(100.0), (10, 2) = 10})

(1)



If you want to use Statistics:-Histogram you must create hypothetic raw data in a correct way.
The reconstruction below produces raw data that have the same histogram as above.

S := [seq((add(op(SyntheticData[k, 1]))/2)$(SyntheticData[k, 2]), k=1..numelems(Classes))]:

Statistics:-Histogram(
  S
  , frequencyscale=absolute
  , binbounds=C
  , axis[1]=[tickmarks=C]
  , view=[-1..101, default]
  , color="Pink"
  , transparency=0.5
  , labels=["Salaires (milliers de \$)", "Nombre d'employés"]
  , labeldirections=[default, vertical]
  , labelfont=[Tahoma, 10]
  , size=[700, 400]
)

 


CASE 3: YOU DO NOT HAVE RAW DATA BUT SYNTHETIC DATA ONLY

... AND YOU DO NOT NECESSARILLY WANT TO USE  Statistics:-Histogram

bars := [
          seq(
            plottools:-rectangle(
              [op(1, SyntheticData[k, 1]), 0]
              , [op(2, SyntheticData[k, 1]), SyntheticData[k, 2]]
              , color="Pink"
              , transparency=0.5
            )
            , k=1..numelems(Classes)
          )
        ]:

plots:-display(
  bars
  , labels=["Salaires (milliers de \$)", "Nombre d'employés"]
  , labeldirections=[default, vertical]
  , labelfont=[Tahoma, 10]
  , size=[700, 400]
)

 

 


 

Download Different_options.mw

It's likely that some typo or hidden character appeared in ode2 (bold red text).
As I never use  the document mode with 2D inputs, the attached file displays two types of fonts: so correct your own file the way I did if this annoys you.

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

with(PDETools):

with(LinearAlgebra):

with(PolynomialIdeals):

with(plots):

declare(f(x), prime = x);

f(x)*`will now be displayed as`*f

 

`derivatives with respect to`*x*`of functions of one variable will now be displayed with '`

 

theta(x)*`will now be displayed as`*theta

 

`derivatives with respect to`*x*`of functions of one variable will now be displayed with '`

(2)

NULL

ode1 := diff(f(x), x, x, x)+(1/2*(1+p*(diff(u[e](x), x))/u[e]))*f(x)*(diff(f(x), x, x))+p*(D(u[e]))(1-(diff(f(x), x))^2)/u[e]-p*((diff(f(x), x))*(diff(f(x), x, x))-(diff(f(x), x))*(diff(f(x), x, x))) = 0;

diff(diff(diff(f(x), x), x), x)+(1/2)*(1+p*(diff(u[e](x), x))/u[e])*f(x)*(diff(diff(f(x), x), x))+p*(D(u[e]))(1-(diff(f(x), x))^2)/u[e] = 0

 

(diff(diff(theta(x), x), x))/Pr+(1/2)*(1+p*(diff(u[e](x), x))/u[e])*f(x)*(diff(theta(x), x))+16*p*a*sigma*L*T[e]^3*(1-theta(x))/(rho*c[p]*u[e]*U[infinity])+U[infinity]^2*p*(diff(f(x), x))*u[e]*(diff(u[e](x), x))/(c[p]*(T[w]-T[e]))-U[infinity]^2*u[e]^2*(diff(diff(f(x), x), x))^2/(c[p]*(T[w]-T[e])) = 0

(3)

 

bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(8) = 1;

f(0) = 0, (D(f))(0) = 0, (D(f))(8) = 1

 

theta(0) = 0, theta(8) = 1

 

U[infinity]*(1-p/L)

(4)

a1 := [a = 1/2, L = 5, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(T[w]-T[e]) = 10, sigma = 5.67*10^(-8), rho = 1, p = .5];

[a = 1/2, L = 5, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(-T[w]+T[e]) = 10, sigma = 0.5670000000e-7, rho = 1, p = .5]

 

[a = 1/2, L = 6, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(-T[w]+T[e]) = 10, sigma = 0.5670000000e-7, rho = 1, p = .5]

 

[a = 1/2, L = 8, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(-T[w]+T[e]) = 10, sigma = 0.5670000000e-7, rho = 1, p = .5]

(5)

b1 := eval(ode2, a1);
b2 := eval(ode2, a2);
b3 := eval(ode2, a3);

1.428571429*(diff(diff(theta(x), x), x))+.5000000000*f(x)*(diff(theta(x), x))+0.3388446214e-1-0.3388446214e-1*theta(x)+0.8067729084e-4*(diff(diff(f(x), x), x))^2 = 0

 

1.428571429*(diff(diff(theta(x), x), x))+.5000000000*f(x)*(diff(theta(x), x))+0.3992205722e-1-0.3992205722e-1*theta(x)+0.8369300576e-4*(diff(diff(f(x), x), x))^2 = 0

 

1.428571429*(diff(diff(theta(x), x), x))+.5000000000*f(x)*(diff(theta(x), x))+0.5204653387e-1-0.5204653387e-1*theta(x)+0.8754046315e-4*(diff(diff(f(x), x), x))^2 = 0

(6)

c1 := dsolve({bcs1, eval(b1, f(x) = x)}, numeric); c1(0);
c2 := dsolve({bcs1, eval(b2, f(x) = x)}, numeric); c2(0);
c3 := dsolve({bcs1, eval(b3, f(x) = x)}, numeric); c3(0);

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(12, {(1) = .0, (2) = .7110816593827365, (3) = 1.4252612944250316, (4) = 2.1450075553641894, (5) = 2.870698607458807, (6) = 3.6016439190865284, (7) = 4.335950433427024, (8) = 5.072381437739985, (9) = 5.809614980858862, (10) = 6.547193597737968, (11) = 7.281757576380688, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = .4938484453872593, (2, 1) = .33592640765199755, (2, 2) = .43888481949800257, (3, 1) = .6120099871592314, (3, 2) = .32852292968068963, (4, 1) = .8036666351004551, (4, 2) = .20577920071677938, (5, 1) = .9150410854740033, (5, 2) = .10721924293183477, (6, 1) = .9688640080289379, (6, 2) = 0.4626167885323632e-1, (7, 1) = .990390459943389, (7, 2) = 0.165006856163208e-1, (8, 1) = .9975131515965905, (8, 2) = 0.4861478032432699e-2, (9, 1) = .999462252457013, (9, 2) = 0.11840633668385558e-2, (10, 1) = .9999041167369538, (10, 2) = 0.23851519861008805e-3, (11, 1) = .9999870438515986, (11, 2) = 0.4007921134972896e-4, (12, 1) = 1.0, (12, 2) = 0.58484376487594675e-5}, datatype = float[8], order = C_order); YP := Matrix(12, 2, {(1, 1) = .4938484453872593, (1, 2) = -0.237191234908843e-1, (2, 1) = .43888481949800257, (2, 2) = -.12498027451545507, (3, 1) = .32852292968068963, (3, 2) = -.1730836385803869, (4, 1) = .20577920071677938, (4, 2) = -.1591461343770899, (5, 1) = .10721924293183477, (5, 2) = -.10974309693499718, (6, 1) = 0.4626167885323632e-1, (6, 2) = -0.59054851436050264e-1, (7, 1) = 0.165006856163208e-1, (7, 2) = -0.25269084092256245e-1, (8, 1) = 0.4861478032432699e-2, (8, 2) = -0.8689730687886674e-2, (9, 1) = 0.11840633668385558e-2, (9, 2) = -0.24203881956517473e-2, (10, 1) = 0.23851519861008805e-3, (10, 2) = -0.5488360802490405e-3, (11, 1) = 0.4007921134972896e-4, (11, 2) = -0.10245379376869548e-3, (12, 1) = 0.58484376487594675e-5, (12, 2) = -0.16375625411613887e-4}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(12, {(1) = .0, (2) = .7110816593827365, (3) = 1.4252612944250316, (4) = 2.1450075553641894, (5) = 2.870698607458807, (6) = 3.6016439190865284, (7) = 4.335950433427024, (8) = 5.072381437739985, (9) = 5.809614980858862, (10) = 6.547193597737968, (11) = 7.281757576380688, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = -0.45502255784570956e-7, (2, 1) = 0.9444480543660673e-8, (2, 2) = -0.32701855724739926e-7, (3, 1) = -0.13111776668557004e-6, (3, 2) = 0.7051671015686677e-7, (4, 1) = -0.10085891808028474e-6, (4, 2) = 0.12193667249947745e-6, (5, 1) = 0.13308574074702087e-6, (5, 2) = -0.7660319092840244e-7, (6, 1) = 0.16436590818059128e-6, (6, 2) = -0.18881505625378492e-6, (7, 1) = -0.20116971988762455e-7, (7, 2) = 0.7685426832740135e-9, (8, 1) = -0.9569009446204095e-7, (8, 2) = 0.13118291896195288e-6, (9, 1) = -0.31100640223178834e-7, (9, 2) = 0.4639384647499807e-7, (10, 1) = 0.1584209066120383e-7, (10, 2) = -0.35974084912581454e-7, (11, 1) = 0.11642475203736646e-7, (11, 2) = -0.2764855504318049e-7, (12, 1) = .0, (12, 2) = -0.1925248521887889e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8881505625378492e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 12, [theta(x), diff(theta(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[theta(x), diff(theta(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8881505625378492e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 12, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, theta(x), diff(theta(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[theta(x), diff(theta(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error  end try end proc

 

[x = 0., theta(x) = HFloat(0.0), diff(theta(x), x) = HFloat(0.4938484453872593)]

 

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(12, {(1) = .0, (2) = .7109194006521007, (3) = 1.4249868973749966, (4) = 2.1446957738695653, (5) = 2.8704181759615643, (6) = 3.6014311331509616, (7) = 4.335815327649273, (8) = 5.07232139321572, (9) = 5.809624985093265, (10) = 6.54727068912088, (11) = 7.281870021183384, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = .49766174438535715, (2, 1) = .3375757993634647, (2, 2) = .44007354626176926, (3, 1) = .6138389834751503, (3, 2) = .32815856063407745, (4, 1) = .8049921933517767, (4, 2) = .2049224933509404, (5, 1) = .9157759787447678, (5, 2) = .1064979573437486, (6, 1) = .9691886486674092, (6, 2) = 0.45847467449322866e-1, (7, 1) = .990506554830634, (7, 2) = 0.16320440331941466e-1, (8, 1) = .9975470109277185, (8, 2) = 0.479979457345354e-2, (9, 1) = .9994703351484253, (9, 2) = 0.1167147574584417e-2, (10, 1) = .9999056808640605, (10, 2) = 0.23475927713373983e-3, (11, 1) = .9999872681820644, (11, 2) = 0.3939868164434801e-4, (12, 1) = 1.0, (12, 2) = 0.5746972025352357e-5}, datatype = float[8], order = C_order); YP := Matrix(12, 2, {(1, 1) = .49766174438535715, (1, 2) = -0.279454400456164e-1, (2, 1) = .44007354626176926, (2, 2) = -.12801162336374522, (3, 1) = .32815856063407745, (3, 2) = -.1744590166938988, (4, 1) = .2049224933509404, (4, 2) = -.15927332083416956, (5, 1) = .1064979573437486, (5, 2) = -.10934646266608739, (6, 1) = 0.45847467449322866e-1, (6, 2) = -0.58651810580892944e-1, (7, 1) = 0.16320440331941466e-1, (7, 2) = -0.25032143866202835e-1, (8, 1) = 0.479979457345354e-2, (8, 2) = -0.8589685100784845e-2, (9, 1) = 0.1167147574584417e-2, (9, 2) = -0.23880431153508012e-2, (10, 1) = 0.23475927713373983e-3, (10, 2) = -0.5405971765590906e-3, (11, 1) = 0.3939868164434801e-4, (11, 2) = -0.10076942378371007e-3, (12, 1) = 0.5746972025352357e-5, (12, 2) = -0.16091521666157516e-4}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(12, {(1) = .0, (2) = .7109194006521007, (3) = 1.4249868973749966, (4) = 2.1446957738695653, (5) = 2.8704181759615643, (6) = 3.6014311331509616, (7) = 4.335815327649273, (8) = 5.07232139321572, (9) = 5.809624985093265, (10) = 6.54727068912088, (11) = 7.281870021183384, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = -0.4573949665850558e-7, (2, 1) = 0.6242874960547203e-8, (2, 2) = -0.3057749989248847e-7, (3, 1) = -0.12966272475504071e-6, (3, 2) = 0.7217465704512478e-7, (4, 1) = -0.940390491554525e-7, (4, 2) = 0.1178112154625901e-6, (5, 1) = 0.1357038240627767e-6, (5, 2) = -0.8040712939043494e-7, (6, 1) = 0.1607319279832156e-6, (6, 2) = -0.18568759788316905e-6, (7, 1) = -0.2282862282805592e-7, (7, 2) = 0.4786547870156748e-8, (8, 1) = -0.94843869899127e-7, (8, 2) = 0.13055152558459717e-6, (9, 1) = -0.2970037963750006e-7, (9, 2) = 0.4434338125388376e-7, (10, 1) = 0.16200082159941674e-7, (10, 2) = -0.36390277880281064e-7, (11, 1) = 0.11561105767115588e-7, (11, 2) = -0.27217683490459804e-7, (12, 1) = .0, (12, 2) = -0.16987712182963502e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8568759788316905e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 12, [theta(x), diff(theta(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[theta(x), diff(theta(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8568759788316905e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 12, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, theta(x), diff(theta(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[theta(x), diff(theta(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error  end try end proc

 

[x = 0., theta(x) = HFloat(0.0), diff(theta(x), x) = HFloat(0.49766174438535715)]

 

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(12, {(1) = .0, (2) = .7105918488760237, (3) = 1.4244328611291581, (4) = 2.1440659423124293, (5) = 2.8698512980998347, (6) = 3.6010004098896413, (7) = 4.335541208202285, (8) = 5.072198845420454, (9) = 5.809644194733494, (10) = 6.5474258600340365, (11) = 7.282096976039094, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = .5052552537006175, (2, 1) = .3408462755204489, (2, 2) = .442422215142914, (3, 1) = .6174560224034431, (3, 2) = .32742681953867314, (4, 1) = .8076074937875736, (4, 2) = .2032242438976165, (5, 1) = .9172227818187347, (5, 2) = .10507373726470326, (6, 1) = .9698264569016594, (6, 2) = 0.45031941250935854e-1, (7, 1) = .9907341917436213, (7, 2) = 0.15966437963481905e-1, (8, 1) = .9976132778134335, (8, 2) = 0.4678915550247515e-2, (9, 1) = .9994861265216441, (9, 2) = 0.11340645752511988e-2, (10, 1) = .9999087319608427, (10, 2) = 0.22742703929550063e-3, (11, 1) = .9999877052438856, (11, 2) = 0.38072228239622816e-4, (12, 1) = 1.0, (12, 2) = 0.5549311882673039e-5}, datatype = float[8], order = C_order); YP := Matrix(12, 2, {(1, 1) = .5052552537006175, (1, 2) = -0.364325736980702e-1, (2, 1) = .442422215142914, (2, 2) = -.13404823355722884, (3, 1) = .32742681953867314, (3, 2) = -.17717619408563023, (4, 1) = .2032242438976165, (4, 2) = -.15951351711336537, (5, 1) = .10507373726470326, (5, 2) = -.10855688752010623, (6, 1) = 0.45031941250935854e-1, (6, 2) = -0.5785531343159658e-1, (7, 1) = 0.15966437963481905e-1, (7, 2) = -0.2456567964351268e-1, (8, 1) = 0.4678915550247515e-2, (8, 2) = -0.8393290947591843e-2, (9, 1) = 0.11340645752511988e-2, (9, 2) = -0.2324700819301251e-2, (10, 1) = 0.22742703929550063e-3, (10, 2) = -0.5244967168305387e-3, (11, 1) = 0.38072228239622816e-4, (11, 2) = -0.9748390992631722e-4, (12, 1) = 0.5549311882673039e-5, (12, 2) = -0.15538073266825603e-4}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(12, {(1) = .0, (2) = .7105918488760237, (3) = 1.4244328611291581, (4) = 2.1440659423124293, (5) = 2.8698512980998347, (6) = 3.6010004098896413, (7) = 4.335541208202285, (8) = 5.072198845420454, (9) = 5.809644194733494, (10) = 6.5474258600340365, (11) = 7.282096976039094, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = -0.4621216507042107e-7, (2, 1) = 0.18871763805575788e-9, (2, 2) = -0.26440976380442328e-7, (3, 1) = -0.12642870795085007e-6, (3, 2) = 0.7507398457807075e-7, (4, 1) = -0.8067853301949186e-7, (4, 2) = 0.10949148171767816e-6, (5, 1) = 0.1404802397309745e-6, (5, 2) = -0.8760075826073426e-7, (6, 1) = 0.15342379606531546e-6, (6, 2) = -0.1792563287988756e-6, (7, 1) = -0.28042403695503607e-7, (7, 2) = 0.12597322677227419e-7, (8, 1) = -0.9307097617997284e-7, (8, 2) = 0.12914164433307568e-6, (9, 1) = -0.26946194672378708e-7, (9, 2) = 0.4029811073645486e-7, (10, 1) = 0.1687782559135978e-7, (10, 2) = -0.37161040827639946e-7, (11, 1) = 0.11392020809816244e-7, (11, 2) = -0.2635259175010777e-7, (12, 1) = .0, (12, 2) = -0.12553181464882845e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.792563287988756e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 12, [theta(x), diff(theta(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[theta(x), diff(theta(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.792563287988756e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 12, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, theta(x), diff(theta(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[theta(x), diff(theta(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error  end try end proc

 

[x = 0., theta(x) = HFloat(0.0), diff(theta(x), x) = HFloat(0.5052552537006175)]

(7)

:

p1 := odeplot(c1, [x, theta(x)], 0 .. 8, colour = green):

display(p1, p2, p3);

 

NULL

Download PGVTN_mmcdara.mw

PS 1: using PDETools:-declare doesn't ease the task of finding where the error is and you must be doubly careful using it.
PS 2: be consistent and declare u[e] too if you consider it as a function of x (which u[e]' seems to mean, even if u[e] is later defined as a constant function).


As explained in the attached file your code does not work with my Maple 2015?
So I copy-pasted your expression of eq5 and used it as the starting point to get eq6.

Your code is in grey, mine on bold brown; the successive transformations of F (the integral term in your eq6) enableobtaoining a standard form from which it clearly appears that the orthogonality condition of Hermite polynomials (Maple uses their physical form) can be used to simplify the expression of F:

Final remark: as _C2 intervenes through its square, there are two real and opposite values of _C2 such that F(n)=1.

restart;

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)


Here is your file NOT EXECUTED in Maple 2015 for PDEtools:-dchange(...) generates an error with Maple 2015

with(PDEtools):
#with(Maplets[Elements]):

eq0:= E1 = (n+1/2)*_h; param := m1 =_h*sigma^2;

E1 = (n+1/2)*_h

 

m1 = _h*sigma^2

(2)

os := expand(isolate((-_h^2/(2*m1))*diff(psi(x),x,x)+(m1*x^2/2)*psi(x)-E1*psi(x)=0,diff(psi(x),x,x)));

diff(diff(psi(x), x), x) = m1^2*x^2*psi(x)/_h^2-2*m1*E1*psi(x)/_h^2

(3)

eq1 := PDEtools:-dchange([x=sqrt(_h/m1)*zeta,psi(x)=f(zeta)*exp(-zeta^2/2)],os,[zeta,f(zeta)],params={m1,_h,E1});
eq2 := expand(isolate(eq1,diff(f(zeta),zeta$2)));
eq3 := simplify(eval(remove(has,convert(rhs(dsolve(eq2,f(zeta))),Hermite),KummerM),eq0)*exp(-zeta^2/2),radical)       assuming zeta::positive;
eq4 := subs(param,simplify(subs(zeta=sqrt(_h/m1)*x,eq3))) ;
eq5 := unapply(eq4,n,x,sigma);
eq6 := solve(int(eq5(n,x,sigma)^2, x=-infinity...infinity)=1,_C2,useassumptions) assuming _C2 :: real, sigma::positive,n::posint;
 

((diff(diff(f(zeta), zeta), zeta))*exp(-(1/2)*zeta^2)-2*(diff(f(zeta), zeta))*zeta*exp(-(1/2)*zeta^2)-f(zeta)*exp(-(1/2)*zeta^2)+f(zeta)*zeta^2*exp(-(1/2)*zeta^2))*m1/_h = m1*zeta^2*f(zeta)*exp(-(1/2)*zeta^2)/_h-2*m1*E1*f(zeta)*exp(-(1/2)*zeta^2)/_h^2

 

diff(diff(f(zeta), zeta), zeta) = -2*E1*f(zeta)/_h+2*(diff(f(zeta), zeta))*zeta+f(zeta)

 

_C2*HermiteH(n, zeta)*exp(-(1/2)*zeta^2)/2^n

 

_C2*HermiteH(n, (1/sigma^2)^(1/2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2)

 

proc (n, x, sigma) options operator, arrow; _C2*HermiteH(n, (1/sigma^2)^(1/2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2) end proc

 

(4)


From now on: I copy-paste the expression of "your" (that I cannot obtain with Maple 2015)

eq5 := proc (n, x, sigma) options operator, arrow; _C2*HermiteH(n, sqrt(1/sigma^2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2) end proc

proc (n, x, sigma) options operator, arrow; _C2*HermiteH(n, sqrt(1/sigma^2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2) end proc

(5)


Here is the integral contained in your relation "eq6"

After some elementary algebra I get the last relation above

F := Int(eq5(n,x,sigma)^2, x=-infinity...infinity) assuming n::posint;
F := combine(F, exp);
F := eval(F, sqrt(1/sigma^2) = 1/sigma);
F := simplify(IntegrationTools:-Expand(F));
F := IntegrationTools:-Change(F, x/sigma=z, z) assuming sigma > 0;
F := simplify(IntegrationTools:-Expand(F));

Int(_C2^2*HermiteH(n, (1/sigma^2)^(1/2)*x)^2*(2^(-n))^2*(exp(-(1/2)*x^2/sigma^2))^2, x = -infinity .. infinity)

 

Int(_C2^2*HermiteH(n, (1/sigma^2)^(1/2)*x)^2*(2^(-n))^2*exp(-x^2/sigma^2), x = -infinity .. infinity)

 

Int(_C2^2*HermiteH(n, x/sigma)^2*(2^(-n))^2*exp(-x^2/sigma^2), x = -infinity .. infinity)

 

_C2^2*(Int(HermiteH(n, x/sigma)^2*exp(-x^2/sigma^2), x = -infinity .. infinity))*4^(-n)

 

_C2^2*(Int(HermiteH(n, z)^2*exp(-z^2)*sigma, z = -infinity .. infinity))*4^(-n)

 

_C2^2*sigma*(Int(HermiteH(n, z)^2*exp(-z^2), z = -infinity .. infinity))*4^(-n)

(6)


From orthogonality of Hermite polynomials (Physical form used in Maple)
(see Wiki)

select(has, indets(F, function), Int)[];
OrthogonalityCondition := % = (2^n*n!*sqrt(Pi))

Int(HermiteH(n, z)^2*exp(-z^2), z = -infinity .. infinity)

 

Int(HermiteH(n, z)^2*exp(-z^2), z = -infinity .. infinity) = 2^n*factorial(n)*Pi^(1/2)

(7)

F := unapply(eval(F, OrthogonalityCondition), n)

proc (n) options operator, arrow; _C2^2*4^(-n)*sigma*2^n*factorial(n)*Pi^(1/2) end proc

(8)


Thus the expression of _C2 for any n

eq6 := unapply( [solve(F(n)=1, _C2)], n)

proc (n) options operator, arrow; [1/(factorial(n)*Pi^(1/2)*2^(-n)*sigma)^(1/2), -1/(factorial(n)*Pi^(1/2)*2^(-n)*sigma)^(1/2)] end proc

(9)

 


 

Download SolnforOS_dsolve_mmcdara.mw



The two derivatives give exactly the same thing a.

Nevertheless here is a way to do what you want derivatives.mw


I_would_do_that.mw

I declared prnt and clr as optional parameters (first one is boolean initialized to true, second one is string initialized to blue), I introduced a control of the color you pass, and I simplified your code a little bit.
Here are the modifications (complete code in the attached file)

spread2:=proc(
               ....
               {clr::string:="Blue"}, 
               {prnt::boolean:=true}
             )
           ....
           local col := substring(StringTools:-Capitalize(clr), 1..1);
           local COL := `if`(col="R", "Red", `if`(col="G", "Green", "Blue"))  
           ....
           if `not`(member(col, {"R", "G", "B"})) then
             error cat("allowed colors are red, green or blue, received ", clr)
           end if:

           if prnt then
             print(cat("Spread 2 [x,y] Points/Vectors wrt origin ",  COL));
           end if:
                 
           if col = "B" then
               return ....

           elif col="G" then
               return ....

           elif col="R" then
               return ....
          end if;
          end proc:



More controls could be necessary for a safer code (what if the denominators in the green and red cases are null?)

restart:

expr := h(arctan(y/x))+arctan(a/b)^2

h(arctan(y/x))+arctan(a/b)^2

(1)

eval(expr, arctan = (u -> arctan((numer, denom)(u))))

h(arctan(y, x))+arctan(a, b)^2

(2)

 

Download arctan.mw

Unfortunately I gave only Maple 2015 and WeibullPlot is not avaliable.

Look to this file and tell me if the plot I draw corresponds or not to what you would like to get?
I'm afraid I did not corerctly understood what is the issue you refer to. 
Here axis 1 is given 4 subticks and axis 2  3 subricks, you can see that the plot is correct.

restart

with(Statistics):

XX := RandomVariable(Weibull(1, .6)):

AA := Sample(XX, 100):

WeibullPlot(AA, reference = false, style = line, gridlines = true, size = [800, 500], axis = [tickmarks = [default, subticks = 4], color = darkgreen])

WeibullPlot(Vector[row](%id = 18446744078086407582), reference = false, style = line, gridlines = true, size = [800, 500], axis = [tickmarks = [default, subticks = 4], color = darkgreen])

(1)

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(2)

F  := unapply(CDF(Weibull(1, 0.6), z), z):

AA := sort(AA):
N  := numelems(AA):

interface(warnlevel=0):
WB := map(n -> [AA[n], solve(F(z)=n/N)], [$1..N-1]):

plots:-display(
  plot(
    WB
    , color="NavyBlue"
    , gridlines=true
    , axis[1]=[mode=log, color="DarkGreen", tickmarks=[subticks=3]]
    , axis[2]=[mode=log, color="DarkGreen", tickmarks=[subticks=4]]
    , size=[700, 500]
    , axes=boxed
  )
  , plot(Mean(Weibull(1, 0.6)), z=(min..max)(AA), color=red, linestyle=3)
  , plot([[Mean(AA), WB[1][2]], [Mean(AA), WB[-1][2]]], z=(min..max)(AA), color=red, linestyle=3)
)

 

Sample_data    := log[10]~(map2(op, 1, WB)):
Reference_data := log[10]~(map2(op, 2, WB)):
FIT := LinearFit([1, x], Reference_data, Sample_data, [x], output=solutionmodule):
print~(FIT:-Results()):

"residualmeansquare" = 0.483263502868184e-2

 

"residualsumofsquares" = .468765597782138

 

"residualstandarddeviation" = 0.695171563621660e-1

 

"degreesoffreedom" = 97

 

"parametervalues" = (Vector(2, {(1) = 0.566005759094454e-1, (2) = .996156636676824}))

 

"parametervector" = (Vector(2, {(1) = 0.566005759094454e-1, (2) = .996156636676824}))

 

"leastsquaresfunction" = 0.566005759094454e-1+.996156636676824*x

 

Typesetting:-mrow(Typesetting:-ms("standarderrors"), Typesetting:-mo("=", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.2777778em", rspace = "0.2777778em"), Typesetting:-maction(Typesetting:-mfenced(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn("0.00770301199087664", mathvariant = "normal"), rowalign = "", columnalign = "", groupalign = "", rowspan = "1", columnspan = "1"), Typesetting:-mtd(Typesetting:-mn("0.00800405156247162", mathvariant = "normal"), rowalign = "", columnalign = "", groupalign = "", rowspan = "1", columnspan = "1"), rowalign = "", columnalign = "", groupalign = ""), align = "axis", rowalign = "baseline", columnalign = "center", groupalign = "{left}", alignmentscope = "true", columnwidth = "auto", width = "auto", rowspacing = "1.0ex", columnspacing = "0.8em", rowlines = "none", columnlines = "none", frame = "none", framespacing = "0.4em 0.5ex", equalrows = "false", equalcolumns = "false", displaystyle = "false", side = "right", minlabelspacing = "0.8em"), mathvariant = "normal", Typesetting:-msemantics = "RowVector", open = "[", close = "]", Typesetting:-msemantics = "RowVector"), actiontype = "rtableaddress", rtableid = "18446744078229114262"))

 

"confidenceintervals" = (Vector(2, {(1) = 0.413122294090191e-1 .. 0.718889224098718e-1, (2) = .980270809958746 .. 1.01204246339490}))

 

"residuals" = (Vector(4, {(1) = ` 1 .. 99 `*Vector[row], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}))

 

"leverages" = (Vector(4, {(1) = ` 1 .. 99 `*Vector[row], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}))

 

"variancecovariancematrix" = (Matrix(2, 2, {(1, 1) = 0.593363937315893e-4, (1, 2) = 0.259631230458869e-4, (2, 1) = 0.259631230458869e-4, (2, 2) = 0.640648414147044e-4}))

 

"internallystandardizedresiduals" = (Vector(4, {(1) = ` 1 .. 99 `*Vector[column], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}))

 

"externallystandardizedresiduals" = (Vector(4, {(1) = ` 1 .. 99 `*Vector[column], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}))

 

"CookDstatistic" = (Vector(4, {(1) = ` 1 .. 99 `*Vector[column], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order}))

 

"AtkinsonTstatistic" = Vector[column](%id = 18446744078229089086)

(3)

p := plots:-display(
  plot([seq([Reference_data[n], Sample_data[n]], n=1..numelems(Sample_data))], color="NavyBlue")
  , plot(
      FIT:-Results("leastsquaresfunction"), x=(min..max)(Reference_data), color="NavyBlue", linestyle=3
    )
  , gridlines=true
):

p := plottools:-transform((x, y) -> [y, x])(p):

plots:-display(
  p
  , axis[1]=[tickmarks=[seq(k=10.0^k, k=-3..1)], color="DarkGreen"]
  , axis[2]=[tickmarks=[seq(k=10.0^k, k=-3..1)], color="DarkGreen"]
  , view=[(min..max)(Sample_data), (min..max)(Reference_data)]
  , size=[700, 500]
  , axes=boxed
)

 

 

Download WeibullPlot_2015.mw

but an interval.

Just do this to see what happens for a given y as x goes to 0.

limit(f(x, y), x=0)
                            /1\             /1\
                   -1 + sin| - | .. 1 + sin| - |
                            \y/             \y/

Maybe doing this will help you understand what happens

limit(f(1/x, 1/y), x=+infinity)
                   -1 + sin(y) .. 1 + sin(y)
limit(sin(x), x=+infinity)
                            -1 .. 1

You can also compute a directional limit 

limit(f(x, lambda*x), x=0)
                            -2 .. 2

In the attached file you will find an exploration of f(x,y) as you get closer and closer to the point (0, 0) (the value h of the slider means the neighborhood of (0, 0) has size 10-h). 
limit.mw


A quick workaround is

plots:-implicitplot(g=1e-8, x = 0 .. 80, y = 0 .. 60, scaling = constrained, grid=[400, 400])


Nevertheless do not gorget that g=0 on the domain  33 <= x <= 58,  0, <= y <= 49.99226668
gnull.mw

(works only is the pattern to search is numeric, if not adaptations have to be done)

restart

#M3 := Matrix(4, 3, [34, 67, 1, 35, 80, 1, 45, 78, 2, 56, 99, 2]):

M3 := LinearAlgebra:-RandomMatrix(10, 5, generator=0..10)

M3 := Matrix(10, 5, {(1, 1) = 0, (1, 2) = 2, (1, 3) = 6, (1, 4) = 5, (1, 5) = 4, (2, 1) = 6, (2, 2) = 10, (2, 3) = 10, (2, 4) = 0, (2, 5) = 3, (3, 1) = 10, (3, 2) = 5, (3, 3) = 5, (3, 4) = 1, (3, 5) = 6, (4, 1) = 2, (4, 2) = 6, (4, 3) = 1, (4, 4) = 0, (4, 5) = 9, (5, 1) = 8, (5, 2) = 5, (5, 3) = 3, (5, 4) = 3, (5, 5) = 5, (6, 1) = 3, (6, 2) = 10, (6, 3) = 4, (6, 4) = 9, (6, 5) = 4, (7, 1) = 10, (7, 2) = 1, (7, 3) = 6, (7, 4) = 8, (7, 5) = 7, (8, 1) = 4, (8, 2) = 2, (8, 3) = 5, (8, 4) = 3, (8, 5) = 8, (9, 1) = 2, (9, 2) = 0, (9, 3) = 3, (9, 4) = 1, (9, 5) = 3, (10, 1) = 7, (10, 2) = 6, (10, 3) = 6, (10, 4) = 2, (10, 5) = 4})

(1)

LookUp := proc(M::Matrix, cr, P, x)
  local S:
  if cr='column' then
    S := select((n -> is(M[n, P]=x)), [$1..numelems(M[..,1])]);
    if S <> [] then
      return M[S, ..];
    else
      error cat("Column ", P, " does not contain ", convert(x, string))
    end if:
  elif cr='row' then
    S := select((n -> is(M[P, n]=x)), [$1..numelems(M[1, ..])]);
    if S <> [] then
      return M[.., S];
    else
      error cat("Row ", P, " does not contain ", convert(x, string))
    end if:
  end if:
end proc:

# Submatrix made of rows whose column "P" contains value "x"

LookUp(M3, 'column', 5, 4);
LookUp(M3, 'column', 2, 1);

Matrix([[0, 2, 6, 5, 4], [3, 10, 4, 9, 4], [7, 6, 6, 2, 4]])

 

Matrix([[10, 1, 6, 8, 7]])

(2)

# Submatrix made of columns whose row "P" contains value "x"


LookUp(M3, 'row', 2, 10);
LookUp(M3, 'row', 1, 4);

Matrix([[2, 6], [10, 10], [5, 5], [6, 1], [5, 3], [10, 4], [1, 6], [2, 5], [0, 3], [6, 6]])

 

Matrix([[4], [3], [6], [9], [5], [4], [7], [8], [3], [4]])

(3)

type(LookUp(M3, 'row', 2, 31), numeric)

Error, (in LookUp) Row 2 does not contain 31

 

 

Download LookUp.mw

Note that I did not pay attention to performance and that the LookUp function may be inefficient for large matrices.

In case you would like to remove the columns/rows which contain the pattern to search, here is variant of the previous file LookUp_variant.mw

... here is a solution  I wonder what you mean by "The result will be a list of 56 models with 5 monomials."?)

It is efficient while the lengths of the two lists are not much larger:

restart:
with(combinat):
L1 := [x^2*y*alpha[1, 11], x*z^2*alpha[2, 15], y^2*z*alpha[3, 17] + y*z*alpha[3, 8]]:
         
L2 := [alpha[i, 0], alpha[i, 1]*x, alpha[i, 2]*y, alpha[i, 3]*z, alpha[i, 4]*x^2, alpha[i, 5]*y*x, alpha[i, 6]*z*x, alpha[i, 7]*y^2, alpha[i, 8]*z*y, alpha[i, 9]*z^2, alpha[i, 10]*x^3, alpha[i, 11]*y*x^2, alpha[i, 12]*z*x^2, alpha[i, 13]*y^2*x, alpha[i, 14]*z*y*x, alpha[i, 15]*z^2*x, alpha[i, 16]*y^3, alpha[i, 17]*z*y^2, alpha[i, 18]*z^2*y, alpha[i, 19]*z^3]:

#-----------------------------------------------------------------

# Build the cartesian product of the two lists

T:=cartprod([L1, L2]):
L1L2 := NULL:
while not T[finished] do 
  L1L2 := L1L2, T[nextvalue]() 
end do:
L1L2 := [L1L2]:

# Keep only the elements which do not have repeated alphas with the same second index.

map(u -> map2(op, 2, [select(type, indets(u), `indexed`)[]]), L1L2):
Keep := zip((u, v) -> if `not`(member(v[-1], {v[1..-2][]})) then add(u) end if, L1L2, %):

print~ (Keep): # to get Keep displayed one element per line

numelems(Keep)
                               56

A_solution.mw

If your matrix has the structure you give, that is +/- a given quantity at some positions and 0 elsewhere, here is a way=

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

b := Matrix(3, 6, [[-1/2, 0, 1/2, 0, 0, 0], [0, 0, 0, -1/2, 0, 1/2], [0, -1/2, -1/2, 1/2, 1/2, 0]]);

b := Matrix(3, 6, {(1, 1) = -1/2, (1, 2) = 0, (1, 3) = 1/2, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = -1/2, (2, 5) = 0, (2, 6) = 1/2, (3, 1) = 0, (3, 2) = -1/2, (3, 3) = -1/2, (3, 4) = 1/2, (3, 5) = 1/2, (3, 6) = 0})

(2)

opb := op(b)[3];

{(1, 1) = -1/2, (1, 3) = 1/2, (2, 4) = -1/2, (2, 6) = 1/2, (3, 2) = -1/2, (3, 3) = -1/2, (3, 4) = 1/2, (3, 5) = 1/2}

(3)

ropb := rhs~(opb)

{-1/2, 1/2}

(4)

b1 := ropb[1] * ``(b / ropb[1])

b1 := -Typesetting[delayDotProduct](1/2, Matrix(3, 6, {(1, 1) = 1, (1, 2) = 0, (1, 3) = -1, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 1, (2, 5) = 0, (2, 6) = -1, (3, 1) = 0, (3, 2) = 1, (3, 3) = 1, (3, 4) = -1, (3, 5) = -1, (3, 6) = 0}), true)

(5)

b2 := ropb[2] * ``(b / ropb[2])

 

Typesetting[delayDotProduct](1/2, Matrix(3, 6, {(1, 1) = -1, (1, 2) = 0, (1, 3) = 1, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = -1, (2, 5) = 0, (2, 6) = 1, (3, 1) = 0, (3, 2) = -1, (3, 3) = -1, (3, 4) = 1, (3, 5) = 1, (3, 6) = 0}), true)

(6)

 

Download For_instance.mw

For more complex matrices there are some ambiguities about what is the "common factor" to put outside of the matrix: For_instance_2.mw

restart

# list of small greek letters
seq(cat(`#mo("&#`, i, `;")`), i=945..969)



# A procedure to generate a random word of L letters
RandomGreekWord := proc(L)
  local r := rand(945..969);
  cat(
       `#mrow(`, 
       seq(cat(`mo("&#`, r(), `;"),`), i=1..L-1),
       cat(`mo("&#`, r(), `;"))`)
  )
end proc:

# an example
RandomGreekWord(10)

Download Greek.mw

Capital greek letters range from 913 to 937.
So all greek letters can be displayed by replacing

i=945..969

by 

i in GreekIndices

where 

GreekIndices := [$913..937, $945..969]

I didn(t clearly understood if you had only 3 set of parameters (your second question) or 27 (3x3x3).
Nevertheless the two cases are adressed in the attached file.

As you load the Statistics package I guessed you intended to use it somewhere. So you will find in this same file a very short illustration of what you could do.

restart; with(LinearAlgebra); with(Statistics)

NULL

Solving an ODE in maple is straigthforward

a1 := .150

.150

(1)

lambda1 := 284

284

(2)

P__b1 := 284

284

(3)

L__Last := 10

10

(4)

ODE1 := diff(P(x), x) = a1*(P__b1-P(x))/(1+lambda1(P(x)))^(2/3)

diff(P(x), x) = 0.5263157895e-3*(284-P(x))*285^(1/3)

(5)

ic1 := P(0) = 40.42

P(0) = 40.42

(6)

sol1 := dsolve({ODE1, ic1}, numeric)

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

(7)

Z := sol1(L__Last)

[x = 10., P(x) = HFloat(48.71221165797189)]

(8)

First question how can I get only the value of P(x) at LLast?

eval(P(x), sol1(L__Last));

select(has, sol1(L__Last), P(x));

select(has, sol1(L__Last), P(x))[];

HFloat(48.71221165797189)

 

[P(x) = HFloat(48.71221165797189)]

 

P(x) = HFloat(48.71221165797189)

(9)

 

Second question how can I solve this ODE for multiple input and finally obtain a vector with all results?

# my prefered way is this one

restart:


ODE1 := diff(P(x), x) = a1*(P__b1-P(x))/(1+lambda1*(P(x)))^(2/3);

params := convert(indets(ODE1, name) minus {x}, list);

sol1 := dsolve({ODE1, P(0) = 40.42}, numeric, parameters=params);

diff(P(x), x) = a1*(P__b1-P(x))/(1+lambda1*P(x))^(2/3)

 

[P__b1, a1, lambda1]

 

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 := [P__b1 = P__b1, a1 = a1, lambda1 = lambda1]; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 24, [( 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..54, {(1) = 1, (2) = 1, (3) = 0, (4) = 0, (5) = 3, (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}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (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..4, {(1) = 40.42, (2) = Float(undefined), (3) = Float(undefined), (4) = 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..1, {(1) = .1}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..1, {(1, 1) = .0}, datatype = float[8], order = C_order), Array(1..1, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0}, datatype = float[8], order = C_order), Array(1..1, {(1) = 0}, datatype = integer[8]), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .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..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .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..1, {(1) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .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..1, {(1) = .0}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..1, {(1, 1) = .0, (2, 0) = .0, (2, 1) = .0, (3, 0) = .0, (3, 1) = .0, (4, 0) = .0, (4, 1) = .0, (5, 0) = .0, (5, 1) = .0, (6, 0) = .0, (6, 1) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = P(x)]`; if Y[1]*Y[4] < -1 then YP[1] := undefined; return 0 end if; YP[1] := Y[3]*(Y[2]-Y[1])*evalf(1/(Y[1]*Y[4]+1)^(2/3)); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = P(x)]`; if Y[1]*Y[4] < -1 then YP[1] := undefined; return 0 end if; YP[1] := Y[3]*(Y[2]-Y[1])*evalf(1/(Y[1]*Y[4]+1)^(2/3)); 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..4, {(1) = 0., (2) = 40.42, (3) = undefined, (4) = undefined}); _vmap := array( 1 .. 1, [( 1 ) = (1)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `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; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch:  end try else try procname(procname("right")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch:  end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error  end try; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _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, P(x)], (4) = [P__b1 = P__b1, a1 = a1, lambda1 = lambda1]}); _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; _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)

data_Pb1 := [284, 283, 352]:
data_a1  := [0.150, 0.152, 0.145]:
data_lam := [38.57, 50.22, 39.83]:

data := [ data_Pb1, data_a1, data_lam ]

[[284, 283, 352], [.150, .152, .145], [38.57, 50.22, 39.83]]

(11)

# To get the 3x3x3=27 different combinations of the three parameters:

use combinat in
  C := cartprod(data):  
  combs := NULL:
  while not C[finished] do
    combs := combs, C[nextvalue]()
  end do:
end use:

combs := [combs]:

# To get the 27 corrresponding values:

L__Last := 10:

RES := table([]):  # a suggestion
for c in combs do
  sol1(parameters = c):
  RES[c] := eval(P(x), sol1(L__Last));
end do:

eval(RES):
  

# Visualizations.
#
# For instance, to plot P(L__last) as a "function" of P__b1 and a1
# "parameterized" by λ1

# Indices of table RES such that λ1 = 38.57
Indices_1 := select((i -> is(op(3, i)=data_lam[1])), [indices(RES, nolist)]):

# Indices of table RES such that λ1 = 50.22
Indices_2 := select((i -> is(op(3, i)=data_lam[2])), [indices(RES, nolist)]):

# Indices of table RES such that λ1 = 39.83
Indices_3 := select((i -> is(op(3, i)=data_lam[3])), [indices(RES, nolist)]):

# Results for λ1 = 39.83
N_1   := numelems(Indices_1):
pts_1 := Matrix(N_1, 3, (i, j) -> `if`(j=3, RES[Indices_1[i]], Indices_1[i][j])):

# Results for λ1 = 39.83
pts_2 := Matrix(N_1, 3, (i, j) -> `if`(j=3, RES[Indices_2[i]], Indices_2[i][j])):
# Results for λ1 = 39.83
pts_3 := Matrix(N_1, 3, (i, j) -> `if`(j=3, RES[Indices_3[i]], Indices_3[i][j])):

plots:-display(
  Statistics:-ScatterPlot3D(pts_1, symbol=solidcircle, symbolsize=25, color=red),
  Statistics:-ScatterPlot3D(pts_2, symbol=solidcircle, symbolsize=25, color=green),
  Statistics:-ScatterPlot3D(pts_3, symbol=solidcircle, symbolsize=25, color=blue)
);

 

ResIndices := [indices(RES, nolist)]:
ResMatrix  := convert(ResIndices, Matrix):
ResMatrix  := `<|>`(ResMatrix, Vector(numelems(combs), i -> RES[ResIndices[i]]));



fit := unapply(Statistics:-LinearFit([1, u, v, w], ResMatrix, [u, v, w]), [u, v, w]):

fit(params[]);

plots:-display(
  Statistics:-ScatterPlot3D(pts_1, symbol=solidcircle, symbolsize=25, color=red),
  Statistics:-ScatterPlot3D(pts_2, symbol=solidcircle, symbolsize=25, color=green),
  Statistics:-ScatterPlot3D(pts_3, symbol=solidcircle, symbolsize=25, color=blue),

  plot3d(fit(P__b1, a1, data_lam[1]), P__b1=(min..max)(data_Pb1), a1=(min..max)(data_a1), style=wireframe, color=red),
  plot3d(fit(P__b1, a1, data_lam[2]), P__b1=(min..max)(data_Pb1), a1=(min..max)(data_a1), style=wireframe, color=green),

  plot3d(fit(P__b1, a1, data_lam[3]), P__b1=(min..max)(data_Pb1), a1=(min..max)(data_a1), style=wireframe, color=blue)
);

ResMatrix := Vector(4, {(1) = ` 27 x 4 `*Matrix, (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

HFloat(39.10075362545055)+HFloat(0.009896105486091737)*P__b1+HFloat(17.595438201866976)*a1-HFloat(0.03830292313118778)*lambda1

 

 

# In case you have just 3 configurations:

ThreeRes  := convert([data_Pb1, data_a1, data_lam], Matrix)^+:
ThreeConf := convert(ThreeRes, listlist):

ResVector := Vector(numelems(ThreeConf)):
for i from 1 to numelems(ThreeConf) do
  sol1(parameters = ThreeConf[i]):
  ResVector[i] := eval(P(x), sol1(L__Last));
end do:

VectorTitle := Vector[row]([params[], P(L__Last)]):

ThreeRes := `<,>`(VectorTitle, < ThreeRes | ResVector >)

ThreeRes := Matrix(4, 4, {(1, 1) = `#msub(mi("P"),mi("b1"))`, (1, 2) = a1, (1, 3) = lambda1, (1, 4) = P(10), (2, 1) = 284, (2, 2) = .150, (2, 3) = 38.57, (2, 4) = 43.0642521778387, (3, 1) = 283, (3, 2) = .152, (3, 3) = 50.22, (3, 4) = 42.6670938150578, (4, 1) = 352, (4, 2) = .145, (4, 3) = 39.83, (4, 4) = 43.6076477820810})

(12)

``

Download SeveralPossibilities.mw

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