mmcdara

3713 Reputation

17 Badges

5 years, 347 days

MaplePrimes Activity


These are answers submitted by mmcdara

Here is an example 

restart:

f := x -> x^2;
g := x -> x^3-1;
h := x -> abs(f(x)-g(x));

u := x -> ['x', 'f(x)', 'g(x)', 'h(x)'];

X := [seq(k, k = 1 .. 5)]:
M := convert([u(x), eval(u~(X))[]], Matrix);

# or maybe
M := convert([eval(u(x)), eval(u~(X))[]], Matrix);

# M can be saved using 
#    ExcelTools:-Export
#    ExportMatrix
#    ...
#
# For instance

n := cat(currentdir(), /"MyMatrix":
ExportMatrix(n, M);



Create_rable_mmcdara.mw

You could also use series(rhs(sol), x=0, 8) instead of asympt in your second example.

For your first example here is the Maple's solution based on a series expansion of the formal solution.
You will see it differs from MMA's from many points and I wasn't capable to transform th Maple's solution into MMA's one.

series.mw

But I guess you already did this yourself?

For the inverse Laplace transform of a function f(s) might exist, f(s) must vanish to 0 faster than
|s^n * f(s)|  (n strictly positive integer) as s tend to infinity.
This is not the case of d(s)


I don't say it's more efficient than what @tomleslie did, but I do believe it is more elegant and closer yo the formal problem:

 

restart:

with(Statistics):

# Be careful:
# In the conventional parameterization of an exponential distribution
# its parameter (p) has the same dimension as the realizations of the random variable.
# Look at this

E := RandomVariable(Exponential(p)):
PDF(E, t)
 

piecewise(t < 0, 0, exp(-t/p)/p)

(1)

# with this is mind your RV should be defined this way

E__X := RandomVariable(Exponential(1/lambda)):
E__Y := RandomVariable(Exponential(1/mu)):
 

# Instead of explicitely write something like
#      assuming lambda > 0, mu > 0
# it's better (IMO) to use the conditions imposed to lambda and mu
# as parameters of exponential distributions:

att := attributes(E__X)[3]: C__X := att:-Conditions[];
att := attributes(E__Y)[3]: C__Y := att:-Conditions[];

0 < 1/lambda

 

0 < 1/mu

(2)

# Then the expected result simply writes

 Probability(E__Y > E__X)  assuming C__X, C__Y

lambda/(lambda+mu)

(3)

 


 

Download The_probabilistic_way.mw

Obviously 7, 8, 12.. are othees among an infinity
 

restart:
sols := isolve(x^2+1=5*k):
S := op~(2, [sols]);

                 [x = 3 - 5 _Z1, x = 2 - 5 _Z1]

seq(op(eval(S, _Z1=n)), n=-10..0):
sort(rhs~([%]))
[2, 3, 7, 8, 12, 13, 17, 18, 22, 23, 27, 28, 32, 33, 37, 38, 42, 

  43, 47, 48, 52, 53]

 


Look to the attached file
proposal.mw

Like @Mac Dude I'm not sure I truly understand what you want.
Maybe something like this?
Maybe_This.mw

Here is a very simple version (continuity corrections do exist but are not implemented here).
Let E some event whose outcomes are 0 and 1

  • n = sample size
  • k = numer of times E=1 is observed amid these n observations
  • p = the hypothetized probability that E=1
  • alpha = the risk of first kind
  • side = either "L" (left tail), "R", right tail or "LR" two-tailed (symmetric)
     

restart

ExactBinomialTest := proc(n, k, p, alpha, side)
  local B, beta, ev, lb, ub, q:
  uses Statistics:
  B    := RandomVariable(Binomial(n, p)):
  ev   := n*p:
  beta := `if`(side="LR", alpha/2, alpha):
  lb   := Quantile(B, beta):
  ub   := Quantile(B, 1-beta):
  #printf("The expected value of k given P=%a and N=%d is %d\n", p, n, ev):


  if k <= ev then
    q := Probability(B <= k, numeric)
  else
    q := Probability(B >= k, numeric)
  end if:
  if k < lb then
    printf("The null hypothesis P0=%a is rejected\n", p):
  elif k > ub then
    printf("The null hypothesis P0=%a is rejected\n", p):
  else
    printf("The null hypothesis P0=%a is accepted\n", p):
  end if:
  printf("p_value = %1.3f\n", q):
end proc:

ExactBinomialTest(100, 41, 1/2, 0.05, "L");
ExactBinomialTest(100, 41, 1/2, 0.05, "LR"):
ExactBinomialTest(100, 41, 1/2, 0.05, "R")

The null hypothesis P0=1/2 is rejected

p_value = 0.044
The null hypothesis P0=1/2 is accepted
p_value = 0.044
The null hypothesis P0=1/2 is rejected
p_value = 0.044

 

ExactBinomialTest(100, 40, 1/3, 0.10, "L");
ExactBinomialTest(100, 40, 1/3, 0.10, "LR"):
ExactBinomialTest(100, 40, 1/3, 0.10, "R")

The null hypothesis P0=1/3 is rejected

p_value = 0.097
The null hypothesis P0=1/3 is accepted
p_value = 0.097
The null hypothesis P0=1/3 is rejected
p_value = 0.097

 

ExactBinomialTest(10^6, 80, 1e-4, 0.10, "L")

The null hypothesis P0=.1e-3 is rejected

p_value = 0.023

 

 

 

Download ExactBinomialTest.mw

I was writting the attached file as @dharr was posting his own answer.

This is basically the same thing, probably more verbose in my case
 

restart; with(plottools); with(plots)

 

In John Stillwell's book "The Four Pillars of Geometry" he declares the following to be the equations of non-Euclidean "lines" where A,B and C are all non-zero


What Maple code will plot this equation for given values of A, B and C?

Equation := A*z*conjugate(z)+B*(z+conjugate(z))+C = 0:

# assumptions : A::real, B::real, C::real

expand(eval(Equation, z=x+I*y)) assuming x::real, y::real;

A*x^2+A*y^2+2*B*x+C = 0

(1)

expand(%/A);
Student:-Precalculus:-CompleteSquare(%, x);
eq := select(has, lhs(%), {x, y}) = -select(`not`@`has`, lhs(%), {x, y})

x^2+y^2+2*B*x/A+C/A = 0

 

(x+B/A)^2+y^2+C/A-B^2/A^2 = 0

 

(x+B/A)^2+y^2 = -C/A+B^2/A^2

(2)

local gamma:
eq := algsubs(C/A=gamma, algsubs(B/A=beta, eq));

Warning, A new binding for the name `gamma` has been created. The global instance of this name is still accessible using the :- prefix, :-`gamma`.  See ?protect for details.

 

(x+beta)^2+y^2 = beta^2-gamma

(3)

# Case rhs(eq) > 0
#
# eq is the equation of a circle of center (-beta, 0) iif gamma verifies

solve(rhs(eq) >= 0, gamma)

{gamma <= beta^2}

(4)

# Case rhs(eq) = 0
#
# the only couple (x,y) which verifies eq when rhs(eq)=0 is

x=-beta, y=0;  # as x, y and beta were assumed real

x = -beta, y = 0

(5)

# Case rhs(eq) < 0
#
# if rhs(eq) < 0, then eq is of the form (x+beta)^2+y^2 = -r^2 with r > 0

sol := solve(eq, y) assuming rhs(eq) < 0, x::real, y::real;

(-2*beta*x-x^2-gamma)^(1/2), -(-2*beta*x-x^2-gamma)^(1/2)

(6)

# -2*beta*x - x^2 - gamma = -2*beta*x - x^2 - gamma + beta^2 - beta^2
# -2*beta*x - x^2 - gamma = -(2*beta*x + x^2 + beta^2) + (beta^2 - gamma)
# -2*beta*x - x^2 - gamma = -(x + beta)^2 + (beta^2 - gamma)
# -2*beta*x - x^2 - gamma < 0  for each real x
#
# thus each member of sol is a complex number which contradicts the assumption y::real.
#
# Case rhs(eq) < 0 has no (x, y) real solution

pl := proc(a, b, c, h)
  local __beta  := b/a:
  local __gamma := c/a:
  local r := __beta^2 - __gamma;
  if r > 0 then
    implicitplot(eval(eq, [beta=__beta, gamma=__gamma]), x=-h..h, y=-h..h, grid=[4*h^2, 4*h^2], scaling=constrained)
  elif r = 0 then
    plot([[-__beta, 0]], style=point, symbol=circle, symbolsize=15, scaling=constrained)
  else
    WARNING("no real solution");
  end if:
end proc:


pl(1, 2, 1, 5);

 

pl(1, 2, 4, 5);

 

pl(1, 2, 5, 5);

Warning, no real solution

 

 


 

Download NonEuclideanLines_mmcdara.mw

 

Curiously the definitions of a vector basis seem to differ in English and in French.

In English a set F of vectors of some space E is a basis if every element of E can be written as a unique linear combination of elements of F.

In French a set ("family") F which verifies this condition is termed a "generating family".
A family F={u, v} is termed "free" if then a*u+b*v=0 ==> a=b=0.
And F is a basis iif it is a generating family and a free family. 

If the number of elements of F equals the dimension of V, then it's enough to veriffy that the family is a generating one or that it is a free one.
Generally it's easier to verify this second condition (in such a case the family is termed "maximally free"):

restart;
u := <1, 1>:
v := <1, 2>:

# Is {u,v} a free family?
cond := lambda*~u + mu*~v =~ <0, 0>:
solve(convert(cond, list), [lambda, mu])
                     [[lambda = 0, mu = 0]]

Correction is in 1D input mode
 

restart

with(LinearAlgebra); A := Matrix(5, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = tau*`&Pi;u`/mu, (3, 4) = 0, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = `&varpi;`*`&Pi;g`/mu}); B := Matrix(5, 5, {(1, 1) = mu, (1, 2) = 0, (1, 3) = 0, (1, 4) = gamma, (1, 5) = 0, (2, 1) = 0, (2, 2) = mu+omega, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = mu+sigma1, (3, 4) = theta, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = alpha2+gamma+mu, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = sigma1+alpha2, (5, 4) = 0, (5, 5) = theta}); C := A.(1/B); Rank(C); evs := Eigenvalues(C); eig := op(`minus`({entries(evs, nolist)}, {0}))

A := Matrix(5, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = tau*`&Pi;u`/mu, (3, 4) = 0, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = varpi*`&Pi;g`/mu})

 

B := Matrix(5, 5, {(1, 1) = mu, (1, 2) = 0, (1, 3) = 0, (1, 4) = gamma, (1, 5) = 0, (2, 1) = 0, (2, 2) = mu+omega, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = mu+sigma1, (3, 4) = theta, (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = alpha2+gamma+mu, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = sigma1+alpha2, (5, 4) = 0, (5, 5) = theta})

 

C := Matrix(5, 5, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = tau*`&Pi;u`/(mu*(mu+sigma1)), (3, 4) = -tau*`&Pi;u`*theta/(mu*(alpha2+gamma+mu)*(mu+sigma1)), (3, 5) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = -varpi*`&Pi;g`*(sigma1+alpha2)/(mu*(mu+sigma1)*theta), (5, 4) = varpi*`&Pi;g`*(sigma1+alpha2)/(mu*(mu+sigma1)*(alpha2+gamma+mu)), (5, 5) = varpi*`&Pi;g`/(mu*theta)})

 

2

 

evs := Vector(5, {(1) = 0, (2) = 0, (3) = 0, (4) = tau*`&Pi;u`/(mu*(mu+sigma1)), (5) = varpi*`&Pi;g`/(mu*theta)})

 

tau*`&Pi;u`/(mu*(mu+sigma1)), varpi*`&Pi;g`/(mu*theta)

(1)

params := {`&Pi;g` = 2.2, `&Pi;u` = 3.4, alpha2 = .33, mu = .2041, omega = .5, sigma1 = .72, tau = .33, theta = .9, varpi = 0.96e-1};

{`&Pi;g` = 2.2, `&Pi;u` = 3.4, alpha2 = .33, mu = .2041, omega = .5, sigma1 = .72, tau = .33, theta = .9, varpi = 0.96e-1}

(2)

plot(eval(eig, remove(has, params, `&varpi;`)), `&varpi;` = 0 .. .5, filled = true, thickness = 4, color = blue, transparency = .4)

Error, (in plot) unexpected options: [11.97669988*varpi, varpi = 0 .. .5]

 

# what are you trying to plot?

eval(eig, remove(has, params, varpi))

5.948820737, 11.97669988*varpi

(3)

# two ways :
#  1/ plot these 2 quantities

plot([eval(eig, remove(has, params, varpi))], varpi = 0 .. .5, filled = true, thickness = 4, color = [blue, red], transparency = .4)

 


#  2/ plot only the last one

plot(eval(eig, remove(has, params, varpi))[2], varpi = 0 .. .5, filled = true, thickness = 4, color = blue, transparency = .4)

 

 


 

Download graph_mmcdara.mw



I used Maple 2015.2

Main problem; the value of delta is not set, I arbitrarily ghoosed delta:=1:

I changed theta[p] into theta__p (not sure this has any influence but I fell safer to use theta__p which is a variable which has noth1ng to do with theta).
The main point is that your equations 1 and 2 do not contain the functtions theta(eta) nor theta__p(eta).
This mean the inititial system can be spilt into 2 sub-systems:

  • sys__12 made of eq1, eq2 and their bcs,
  • sys__34 made of eq3, eq4 and their bcs.

Thus:

  1. dsolve sys__12 (I use the option output=procedurelist)
     
  2. Define a procedure K wich takes as argument the name eta or a numerical value of eta and returns in this case the evaluations of [f(eta), F(eta), diff(f(eta), eta), diff(F(eta), eta)] for this value of theta.
     
  3. After having instanciate sys__34 with params, replace:
    1. f(eta) by K(eta)[1],
    2. F(eta) by K(eta)[2],
    3. diff(f(eta), eta) by K(eta)[3],
    4. diff(F(eta), eta) by K(eta)[4],
       
  4. Solve the resulting sys__34 system.
    You must use the option known=K to tell Maple that K(eta)[1]...K(eta)[4] are known elsewhere.know

    Without any precaution you should get ab error which advices you tou use the midpoint submethod either midrich or middefer.
    Whatever the submethod I used the system seems so badly conditionned that I wasn't able to reach a reasonable accuracy (I desperately set abserr=1 to obtain a solution). This what I titled my answer "almost done".

restart

with(plots):

eq1 := (2*eta*gamma+1)*(diff(f(eta), `$`(eta, 3)))+2*gamma*(diff(f(eta), `$`(eta, 2)))+f(eta)*(diff(f(eta), `$`(eta, 2)))-(diff(f(eta), eta))^2-(Q+S)*(diff(f(eta), eta))+beta*(diff(F(eta), eta)-(diff(f(eta), eta))) = 0:

eq2 := (diff(F(eta), `$`(eta, 2)))*F(eta)-(diff(F(eta), eta))^2+beta*(diff(f(eta), eta)-(diff(F(eta), eta))) = 0:

eq3 := (2*eta*gamma+1)*(1+Rd)*(diff(theta(eta), `$`(eta, 2))) = -Pr*((diff(theta(eta), eta))*f(eta)-2*(diff(f(eta), eta))*theta(eta))-gamma*(diff(theta(eta), eta))-N*Pr*betat*((theta__p(eta), eta)-theta(eta))-N*Pr*Ec*betat*(diff(F(eta), eta)-(diff(f(eta), eta)))-Pr*delta*theta(eta):

eq4 := 2*(diff(theta__p(eta), eta))*f(eta)-F(eta)*theta__p(eta)+betat*delta*(theta__p(eta)-theta(eta)) = 0:

bcs := f(0) = 0, (D(f))(0) = 1, (D(f))(5) = 0, (D(F))(5) = 0, F(5) = f(5), theta(0) = 1, theta(5) = 0, theta__p(5) = 0:

params := [Rd = .1, beta = .5, Q = .5, S = .5, gamma = .1, Pr = 6.2, N = .5, betat = .5, Ec = .1]:

delta := 1:  # delta not fixed !!!

# Main observation; eq1 and eq2 depend only upon f and F  

indets(eq1, function);
indets(eq2, function);

{F(eta), diff(F(eta), eta), diff(diff(diff(f(eta), eta), eta), eta), diff(diff(f(eta), eta), eta), diff(f(eta), eta), f(eta)}

 

{F(eta), diff(F(eta), eta), diff(diff(F(eta), eta), eta), diff(f(eta), eta), f(eta)}

(1)

sys__12 := eval([eq1, eq2, f(0) = 0, (D(f))(0) = 1, (D(f))(5) = 0, (D(F))(5) = 0, F(5) = f(5)], params):
sol__12 := dsolve(sys__12, numeric, output = procedurelist);

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(17, {(1) = .0, (2) = .2700295406680469, (3) = .5461814747766728, (4) = .8292588209972042, (5) = 1.1234150774760683, (6) = 1.431035339629634, (7) = 1.7496291223773444, (8) = 2.075914881269067, (9) = 2.408230078854766, (10) = 2.742205295476464, (11) = 3.0774349455365573, (12) = 3.4129672931189416, (13) = 3.74869641693143, (14) = 4.084468803510717, (15) = 4.417149968189033, (16) = 4.721816642175481, (17) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(17, 5, {(1, 1) = .44776527313438064, (1, 2) = .33133727618764913, (1, 3) = .0, (1, 4) = 1.0, (1, 5) = -1.5719460272279495, (2, 1) = .5213624951662373, (2, 2) = .2217023909765899, (2, 3) = .22045952998207696, (2, 4) = .6583232859380348, (2, 5) = -1.0037671287227639, (3, 1) = .5718234259032432, (3, 2) = .14894670237473073, (3, 3) = .36912084196073325, (3, 4) = .43487740306649103, (3, 5) = -.6430155320436061, (4, 1) = .6066140717509945, (4, 2) = .10029249480024327, (4, 3) = .469902127251469, (4, 4) = .2879720893049211, (4, 5) = -.4130037744612695, (5, 1) = .6309127884296852, (5, 2) = 0.6725940922492528e-1, (5, 3) = .5391442607170165, (5, 4) = .19006598529443805, (5, 5) = -.26445645144099067, (6, 1) = .6478941036296227, (6, 2) = 0.4475730981457941e-1, (6, 3) = .586800586045253, (6, 4) = .12467762652820233, (6, 5) = -.1683906813056968, (7, 1) = .6595684318193582, (7, 2) = 0.29614870473976596e-1, (7, 3) = .6191402504248181, (7, 4) = 0.8157762164726536e-1, (7, 5) = -.10711429298647333, (8, 1) = .6674672167892126, (8, 2) = 0.1951756867756001e-1, (8, 3) = .6408263407473044, (8, 4) = 0.5344481967149114e-1, (8, 5) = -0.6842360971338068e-1, (9, 1) = .6727561100925025, (9, 2) = 0.12782783201441923e-1, (9, 3) = .6553122229469226, (9, 4) = 0.3508354643667022e-1, (9, 5) = -0.44010428510442734e-1, (10, 1) = .6762273833779221, (10, 2) = 0.8307953068430345e-2, (10, 3) = .6648932981595975, (10, 4) = 0.23142480667606512e-1, (10, 5) = -0.28681021784476287e-1, (11, 1) = .6784756350927372, (11, 2) = 0.5303277034954629e-2, (11, 3) = .6712427107983656, (11, 4) = 0.15279463138533988e-1, (11, 5) = -0.18963855581152253e-1, (12, 1) = .679892179218728, (12, 2) = 0.32714663294575444e-2, (12, 3) = .6754314367159898, (12, 4) = 0.10033737564104571e-1, (12, 5) = -0.12763512056209856e-1, (13, 1) = .6807440254039239, (13, 2) = 0.18925878798273792e-2, (13, 3) = .6781638938053991, (13, 4) = 0.6466085693279025e-2, (13, 5) = -0.878177393716817e-2, (14, 1) = .6812136245064051, (14, 2) = 0.9684873470550855e-3, (14, 3) = .6798934576626122, (14, 4) = 0.3978616235233124e-2, (14, 5) = -0.62238500751895555e-2, (15, 1) = .6814306856085619, (15, 2) = 0.3850239255045409e-3, (15, 3) = .6809058823403, (15, 4) = 0.21970566082543854e-2, (15, 5) = -0.4609984823372287e-2, (16, 1) = .6814973367424689, (16, 2) = 0.8800363841388453e-4, (16, 3) = .6813773865384698, (16, 4) = 0.9452622479871761e-3, (16, 5) = -0.36799735133912046e-2, (17, 1) = .6815055279510432, (17, 2) = .0, (17, 3) = .6815055279510432, (17, 4) = .0, (17, 5) = -0.3161424611013035e-2}, datatype = float[8], order = C_order); YP := Matrix(17, 5, {(1, 1) = .33133727618764913, (1, 2) = -.5014836674199502, (1, 3) = 1.0, (1, 4) = -1.5719460272279495, (1, 5) = 2.648720567351765, (2, 1) = .2217023909765899, (2, 2) = -.3244546719112376, (2, 3) = .6583232859380348, (2, 4) = -1.0037671287227639, (2, 5) = 1.643317863794819, (3, 1) = .14894670237473073, (3, 2) = -.21121945119123267, (3, 3) = .43487740306649103, (3, 4) = -.6430155320436061, (3, 5) = 1.0213465383823663, (4, 1) = .10029249480024327, (4, 2) = -.1381128738033511, (4, 3) = .2879720893049211, (4, 4) = -.4130037744612695, (4, 5) = .635940126024083, (5, 1) = 0.6725940922492528e-1, (5, 2) = -0.9015423517890793e-1, (5, 3) = .19006598529443805, (5, 4) = -.26445645144099067, (5, 5) = .3944415120283217, (6, 1) = 0.4475730981457941e-1, (6, 2) = -0.5858510111811698e-1, (6, 3) = .12467762652820233, (6, 4) = -.1683906813056968, (6, 5) = .24309630227876206, (7, 1) = 0.29614870473976596e-1, (7, 2) = -0.3806175951175528e-1, (7, 3) = 0.8157762164726536e-1, (7, 4) = -.10711429298647333, (7, 5) = .1496049121417643, (8, 1) = 0.1951756867756001e-1, (8, 2) = -0.2484420147202396e-1, (8, 3) = 0.5344481967149114e-1, (8, 4) = -0.6842360971338068e-1, (8, 5) = 0.9242420908739084e-1, (9, 1) = 0.12782783201441923e-1, (9, 2) = -0.1633130031287029e-1, (9, 3) = 0.3508354643667022e-1, (9, 4) = -0.44010428510442734e-1, (9, 5) = 0.57441143040287745e-1, (10, 1) = 0.8307953068430345e-2, (10, 2) = -0.10866524923457799e-1, (10, 3) = 0.23142480667606512e-1, (10, 4) = -0.28681021784476287e-1, (10, 5) = 0.3610169213484274e-1, (11, 1) = 0.5303277034954629e-2, (11, 2) = -0.7310458987672048e-2, (11, 3) = 0.15279463138533988e-1, (11, 4) = -0.18963855581152253e-1, (11, 5) = 0.229176337398394e-1, (12, 1) = 0.32714663294575444e-2, (12, 2) = -0.4957305332224504e-2, (12, 3) = 0.10033737564104571e-1, (12, 4) = -0.12763512056209856e-1, (12, 5) = 0.14673258499433704e-1, (13, 1) = 0.18925878798273792e-2, (13, 2) = -0.3353928837624718e-2, (13, 3) = 0.6466085693279025e-2, (13, 4) = -0.878177393716817e-2, (13, 5) = 0.9433680672377034e-2, (14, 1) = 0.9684873470550855e-3, (14, 2) = -0.2208009972550792e-2, (14, 3) = 0.3978616235233124e-2, (14, 4) = -0.62238500751895555e-2, (14, 5) = 0.6040988893516817e-2, (15, 1) = 0.3850239255045409e-3, (15, 2) = -0.13293620570413733e-2, (15, 3) = 0.21970566082543854e-2, (15, 4) = -0.4609984823372287e-2, (15, 5) = 0.38062804457790473e-2, (16, 1) = 0.8800363841388453e-4, (16, 2) = -0.6289409173556955e-3, (16, 3) = 0.9452622479871761e-3, (16, 4) = -0.36799735133912046e-2, (16, 5) = 0.2375189062515268e-2, (17, 1) = .0, (17, 2) = .0, (17, 3) = .0, (17, 4) = -0.3161424611013035e-2, (17, 5) = 0.13934066354042334e-2}, 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(17, {(1) = .0, (2) = .2700295406680469, (3) = .5461814747766728, (4) = .8292588209972042, (5) = 1.1234150774760683, (6) = 1.431035339629634, (7) = 1.7496291223773444, (8) = 2.075914881269067, (9) = 2.408230078854766, (10) = 2.742205295476464, (11) = 3.0774349455365573, (12) = 3.4129672931189416, (13) = 3.74869641693143, (14) = 4.084468803510717, (15) = 4.417149968189033, (16) = 4.721816642175481, (17) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(17, 5, {(1, 1) = -0.5557170630311661e-9, (1, 2) = -0.8568444294619354e-8, (1, 3) = .0, (1, 4) = .0, (1, 5) = -0.996321158246133e-8, (2, 1) = -0.29775367430238932e-7, (2, 2) = 0.51448201034330047e-7, (2, 3) = -0.9338503564926378e-7, (2, 4) = 0.16119701358192726e-6, (2, 5) = -0.2661709248184363e-6, (3, 1) = -0.2891050343668424e-7, (3, 2) = 0.5344899476525872e-7, (3, 3) = -0.9176042528308786e-7, (3, 4) = 0.15799013994899224e-6, (3, 5) = -0.2499917438861638e-6, (4, 1) = -0.1927807403261658e-7, (4, 2) = 0.38893522217972736e-7, (4, 3) = -0.6314408631390969e-7, (4, 4) = 0.11040192879654927e-6, (4, 5) = -0.168450566115583e-6, (5, 1) = -0.9451108871656315e-8, (5, 2) = 0.23506561485624035e-7, (5, 3) = -0.3433544296514957e-7, (5, 4) = 0.6357224004387666e-7, (5, 5) = -0.9346927926922444e-7, (6, 1) = -0.17952890062909867e-8, (6, 2) = 0.11485578153704378e-7, (6, 3) = -0.12371911215185314e-7, (6, 4) = 0.2856188378954896e-7, (6, 5) = -0.40177124898242385e-7, (7, 1) = 0.33629876485582894e-8, (7, 2) = 0.3438144666793098e-8, (7, 3) = 0.19670495465459734e-8, (7, 4) = 0.6170587583276722e-8, (7, 5) = -0.7855938356527589e-8, (8, 1) = 0.6290794653972546e-8, (8, 2) = -0.10736074960720638e-8, (8, 3) = 0.969080021609951e-8, (8, 4) = -0.569460732839434e-8, (8, 5) = 0.8253813531999202e-8, (9, 1) = 0.7572046363127913e-8, (9, 2) = -0.3010207778743067e-8, (9, 3) = 0.12689878914657082e-7, (9, 4) = -0.1033127198991506e-7, (9, 5) = 0.13987204260214792e-7, (10, 1) = 0.7829327930419819e-8, (10, 2) = -0.3347426371048808e-8, (10, 3) = 0.12839070161904223e-7, (10, 4) = -0.1073328829776e-7, (10, 5) = 0.14076279806383137e-7, (11, 1) = 0.7582703863603602e-8, (11, 2) = -0.2891392769667636e-8, (11, 3) = 0.11656345534185066e-7, (11, 4) = -0.9198752266915943e-8, (11, 5) = 0.11835993527707034e-7, (12, 1) = 0.7166421536865248e-8, (12, 2) = -0.2148435088658628e-8, (12, 3) = 0.1007706561001513e-7, (12, 4) = -0.7061367623274055e-8, (12, 5) = 0.9048774297591804e-8, (13, 1) = 0.67611797018740185e-8, (13, 2) = -0.13976774396261053e-8, (13, 3) = 0.8595435301400648e-8, (13, 4) = -0.4983619985985624e-8, (13, 5) = 0.6526712448285493e-8, (14, 1) = 0.6440901841297378e-8, (14, 2) = -0.7706065262608924e-9, (14, 3) = 0.741967636293521e-8, (14, 4) = -0.3217163270911993e-8, (14, 5) = 0.454060837838052e-8, (15, 1) = 0.6216116456491702e-8, (15, 2) = -0.320565246765447e-9, (15, 3) = 0.6605099271231325e-8, (15, 4) = -0.18085873364690693e-8, (15, 5) = 0.3117357572183856e-8, (16, 1) = 0.6064304046619167e-8, (16, 2) = -0.6701463087834589e-10, (16, 3) = 0.6136443374376152e-8, (16, 4) = -0.7630971430389972e-9, (16, 5) = 0.2232644899187691e-8, (17, 1) = 0.5954749976109993e-8, (17, 2) = .0, (17, 3) = 0.5954749976109993e-8, (17, 4) = .0, (17, 5) = 0.17510675834850087e-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[17] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.661709248184363e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [5, 17, [F(eta), diff(F(eta), eta), f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[17] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[17] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(5, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(17, 5, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(5, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(17, 5, X, Y, outpoint, yout, L, V) end if; [eta = outpoint, seq('[F(eta), diff(F(eta), eta), f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta)]'[i] = yout[i], i = 1 .. 5)] 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[17] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.661709248184363e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [5, 17, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[17] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[17] 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(5, {(1) = .0, (2) = .0, (3) = .0, (4) = .0, (5) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(17, 5, 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(5, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0.}); `dsolve/numeric/hermite`(17, 5, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 5)] end proc, (2) = Array(0..0, {}), (3) = [eta, F(eta), diff(F(eta), eta), f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta)], (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); [eta = res[1], seq('[F(eta), diff(F(eta), eta), f(eta), diff(f(eta), eta), diff(diff(f(eta), eta), eta)]'[i] = res[i+1], i = 1 .. 5)] catch: error  end try end proc

 

[eta = 2., F(eta) = HFloat(0.6659115637409204), diff(F(eta), eta) = HFloat(0.021498918991851815), f(eta) = HFloat(0.6365649907511923), diff(f(eta), eta) = HFloat(0.058915598688308896), diff(diff(f(eta), eta), eta) = HFloat(-0.07584262138836992)]

(2)

K := proc(s)
  local Kf, KF, Kdf, KdF:
  if not type(evalf(s),'numeric') then
    'procname'(s);
  else
    return eval([f(eta), F(eta), diff(f(eta), eta), diff(F(eta), eta)], sol__12(s))
  end if:
end proc:
  
# Check
plot(['K(s)[1]', 'K(s)[2]'], s=1..5, color=[red, blue], legend=['f(eta)', 'F(eta)'], labels=[eta, ""]);
 

 

sys__34 := isolate(eq3, diff(theta(eta), eta, eta)),
           isolate(eq4, diff(theta__p(eta), eta)):
sys__34 := eval([sys__34, theta(0) = 1, theta(5) = 0, theta__p(5) = 0], params);

[diff(diff(theta(eta), eta), eta) = .9090909091*(-6.2*(diff(theta(eta), eta))*f(eta)+12.4*(diff(f(eta), eta))*theta(eta)-.1*(diff(theta(eta), eta))-1.550*theta__p(eta)-4.650*theta(eta)-.1550*(diff(F(eta), eta))+.1550*(diff(f(eta), eta)))/(.2*eta+1), diff(theta__p(eta), eta) = (1/2)*(F(eta)*theta__p(eta)-.5*theta__p(eta)+.5*theta(eta))/f(eta), theta(0) = 1, theta(5) = 0, theta__p(5) = 0]

(3)

sys__34 := eval(sys__34, [f(eta)='K(eta)[1]', F(eta)='K(eta)[2]', diff(f(eta), eta)='K(eta)[3]', diff(F(eta), eta)='K(eta)[4]']);

[diff(diff(theta(eta), eta), eta) = .9090909091*(-6.2*(diff(theta(eta), eta))*K(eta)[1]+12.4*K(eta)[3]*theta(eta)-.1*(diff(theta(eta), eta))-1.550*theta__p(eta)-4.650*theta(eta)-.1550*K(eta)[4]+.1550*K(eta)[3])/(.2*eta+1), diff(theta__p(eta), eta) = (1/2)*(K(eta)[2]*theta__p(eta)-.5*theta__p(eta)+.5*theta(eta))/K(eta)[1], theta(0) = 1, theta(5) = 0, theta__p(5) = 0]

(4)

sol__34 := dsolve(sys__34, numeric, output=procedurelist, known=K, method=bvp[middefer], maxmesh=512, abserr=1e-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(11, {(1) = .0, (2) = .3370404083306422, (3) = .7008767283468323, (4) = 1.1112972869729818, (5) = 1.609997064526372, (6) = 2.189609762048589, (7) = 2.812077057741318, (8) = 3.44210959171944, (9) = 4.073127629221379, (10) = 4.627447082063739, (11) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(11, 3, {(1, 1) = 1.0, (1, 2) = -2.631180930496676, (1, 3) = -.5069651588422761, (2, 1) = .4023510464572303, (2, 2) = -.9152718910087276, (2, 3) = -.10719058799533322, (3, 1) = .1726623490489365, (3, 2) = -.34732166973349327, (3, 3) = -0.3642234782724219e-1, (4, 1) = 0.7245529317024713e-1, (4, 2) = -.14099234754733977, (4, 3) = -0.11955501922042339e-1, (5, 1) = 0.23698556314753913e-1, (5, 2) = -0.5454307897597913e-1, (5, 3) = -0.20139753068670863e-2, (6, 1) = 0.3049194442117786e-2, (6, 2) = -0.16709196762319765e-1, (6, 3) = 0.10235425366078611e-2, (7, 1) = -0.2719856786139026e-2, (7, 2) = -0.18268814161701507e-2, (7, 3) = 0.11518671538119254e-2, (8, 1) = -0.2614649585026932e-2, (8, 2) = 0.2160855283988688e-2, (8, 3) = 0.600072097033077e-3, (9, 1) = -0.12708908932571583e-2, (9, 2) = 0.20981630384967246e-2, (9, 3) = 0.18102423815965176e-3, (10, 1) = -0.3515015852953094e-3, (10, 2) = 0.12190191503928644e-2, (10, 3) = 0.23439573975186445e-4, (11, 1) = .0, (11, 2) = 0.6679696148136223e-3, (11, 3) = .0}, datatype = float[8], order = C_order); YP := Matrix(11, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = -.9152718910087276, (2, 2) = 2.4727536746209924, (2, 3) = .3760823990110259, (3, 1) = -.34732166973349327, (3, 2) = .7893743028047262, (3, 3) = 0.9662576107625089e-1, (4, 1) = -.14099234754733977, (4, 2) = .2662816956052928, (4, 3) = 0.3229410448051668e-1, (5, 1) = -0.5454307897597913e-1, (5, 2) = 0.9768863241923123e-1, (5, 3) = 0.9509064038998313e-2, (6, 1) = -0.16709196762319765e-1, (6, 2) = 0.37412194655996434e-1, (6, 3) = 0.13133552399215123e-2, (7, 1) = -0.18268814161701507e-2, (7, 2) = 0.1163249008646273e-1, (7, 3) = -0.867516162778065e-3, (8, 1) = 0.2160855283988688e-2, (8, 2) = 0.14317512514059364e-2, (8, 3) = -0.8874409868549002e-3, (9, 1) = 0.20981630384967246e-2, (9, 2) = -0.15103520383209646e-2, (9, 3) = -0.44321926460382834e-3, (10, 1) = 0.12190191503928644e-2, (10, 2) = -0.16528238604940286e-2, (10, 3) = -0.12586527735752815e-3, (11, 1) = 0.6679696148136223e-3, (11, 2) = -0.13132690311376787e-2, (11, 3) = .0}, 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(11, {(1) = .0, (2) = .3370404083306422, (3) = .7008767283468323, (4) = 1.1112972869729818, (5) = 1.609997064526372, (6) = 2.189609762048589, (7) = 2.812077057741318, (8) = 3.44210959171944, (9) = 4.073127629221379, (10) = 4.627447082063739, (11) = 5.0}, datatype = float[8], order = C_order); Y := Matrix(11, 3, {(1, 1) = .0, (1, 2) = -0.13944241604693985e-1, (1, 3) = 0.2199320251397907e-2, (2, 1) = 0.3122582429547549e-1, (2, 2) = -.10070764912778984, (2, 3) = -0.10906300241481343e-1, (3, 1) = 0.2396545586162674e-1, (3, 2) = -0.8321640362333171e-1, (3, 3) = 0.3933131209080537e-3, (4, 1) = 0.8640618859978494e-2, (4, 2) = -0.4791621713263419e-1, (4, 3) = 0.3229350888215585e-2, (5, 1) = -0.17331011525628486e-2, (5, 2) = -0.1902105958328692e-1, (5, 3) = 0.32869341286691613e-2, (6, 1) = -0.4847665589688804e-2, (6, 2) = -0.29053410124298754e-2, (6, 3) = 0.21332021422819417e-2, (7, 1) = -0.36149600965167117e-2, (7, 2) = 0.197695795668519e-2, (7, 3) = 0.9832255075634968e-3, (8, 1) = -0.17708414039592315e-2, (8, 2) = 0.2040124287083247e-2, (8, 3) = 0.32939453342591216e-3, (9, 1) = -0.60979066315007e-3, (9, 2) = 0.11477551977751616e-2, (9, 3) = 0.680826223469885e-4, (10, 1) = -0.13922033932913127e-3, (10, 2) = 0.5435456008739585e-3, (10, 3) = 0.5485745412833969e-5, (11, 1) = .0, (11, 2) = 0.29761609702088786e-3, (11, 3) = .0}, 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[11] elif outpoint = "order" then return 2 elif outpoint = "error" then return .100707649127789844 elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [3, 11, [theta(eta), diff(theta(eta), eta), theta__p(eta)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[11] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[11] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(3, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(11, 3, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(3, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(11, 3, X, Y, outpoint, yout, L, V) end if; [eta = outpoint, seq('[theta(eta), diff(theta(eta), eta), theta__p(eta)]'[i] = yout[i], i = 1 .. 3)] 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[11] elif outpoint = "order" then return 2 elif outpoint = "error" then return .100707649127789844 elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [3, 11, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[11] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[11] 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 ) = (true), ( 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(3, {(1) = .0, (2) = .0, (3) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(11, 3, 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 ) = (true), ( 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(3, {(1) = 0., (2) = 0., (3) = 0.}); `dsolve/numeric/hermite`(11, 3, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 3)] end proc, (2) = Array(0..0, {}), (3) = [eta, theta(eta), diff(theta(eta), eta), theta__p(eta)], (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); [eta = res[1], seq('[theta(eta), diff(theta(eta), eta), theta__p(eta)]'[i] = res[i+1], i = 1 .. 3)] catch: error  end try end proc

(5)

display(
  odeplot(sol__34, [eta, theta(eta)], eta=0..5, color=red, legend='theta(eta)'),
  odeplot(sol__34, [eta, theta__p(eta)], eta=0..5, color=blue, legend='theta__p(eta)')
)

 

 


 

Download MapleOde_mmcdara.mw

According to my previous reply (correct expression of bcs and ics).
You 'll see in the attached file that there is no more RootOf.

elctroWav_mmcdara.mw

I'm ,evertheless a little bit surprised by your boundary conditions.
What are R and rb?
It seems that rho varies between 0 and 1 and that R-rb is the inner radius and R the outer radius of a tube within which a wave propagates?
Why don't set dyrectly bcs at rho=R-rb and rho=R?
 

restart
interface(version)
Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 

   21 2015 Build ID 1097895
p := .5:
n := 10:
nseq := 2:
with(Statistics):
randomize():
UseHardwareFloats := false:
x1 := seq(Sample(RandomVariable(BernoulliDistribution(p)), n), i = 1 .. nseq):
x2 := seq(convert(x1[i], list), i = 1 .. nseq);
        [0., 1.0, 0., 0., 1.0, 0., 0., 0., 1.0, 1.0], 

          [1.0, 1.0, 1.0, 1.0, 0., 1.0, 0., 0., 0., 0.]
numboccur(x2[1], {0.})
                               6

An alternative to numboccur:
 

Tally(x2[1])
                       [0. = 6, 1.0 = 4]

 

How to find the answer to question c with MAPLE?
 

restart:
with(IntegrationTools):
J16 := Int(f(x), x=1..6);
J16 = Split(J16, 4);
isolate(%, Int(f(x), x=4..6));
eval(%, [Int(f(x), x=1..4) = 7, Int(f(x), x=1..6) = 4])

                         /6             
                        |               
                        |   f(x) dx = -3
                        |               
                       /4               


Extra: Answer to question f could be obtained this way

Qf := Int(3*f(x)-6*g(x), x=1..6):
Qf = Expand(Qf):
lhs(%) = eval(rhs(%), [Int(f(x), x=1..6) = 4, Int(g(x), x=1..6) = 5]);
                   /6                         
                  |                           
                  |   3 f(x) - 6 g(x) dx = -18
                  |                           
                 /1                           

 

3 4 5 6 7 8 9 Last Page 5 of 30