tomleslie

4487 Reputation

10 Badges

9 years, 98 days

MaplePrimes Activity


These are answers submitted by tomleslie

you have two problems, one syntactic and one logical

Syntatic

In Maple you can run into all sorts of weird problems if you use the same name for indexed and unindexed variables - so writing anything like f(x):=sum(f[i](x),i=1..N) is just asking for trouble. Where you have done this I have capitalised the left-hand-side to produce distict names. (Obviously, I have appropriately corrected all subsequent name references in the code)

Logical

Becuase your equations are coupled, you cannot solve for f[i] i=0..N, then theta[i] i=0..N and so on sequentially. However for any value of 'i' you can simultaneously solve for the [ f[i], theta[i], u[i], w[i] ], since these only depend on previously obtained values, ie [ f[i-1], theta[i-1], u[i-1], w[i-1] ].

The attached now executes withut error, and appears to return all required quantities Be aware that some of these expressions are a bit "lengthy"


 

  restart;
  PDEtools[declare](f(x),theta(x),u(x),w(x), prime=x):
  N:= 4;

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

 

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

 

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

 

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

 

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

 

4

(1)

  F:=      x -> local i;
                add
                ( p^i*f[i](x),
                  i=0..N
                );
  Theta:=  x -> local i;
                add
                ( p^i*theta[i](x),
                  i = 0..N
                );
  U:=      x -> local i;
                add
                ( p^i*u[i](x),
                  i = 0..N
                );
  W:=      x -> local i;
                add
                ( p^i*w[i](x),
                  i = 0 .. N
                );

proc (x) local i; options operator, arrow; add(p^i*f[i](x), i = 0 .. N) end proc

 

proc (x) local i; options operator, arrow; add(p^i*theta[i](x), i = 0 .. N) end proc

 

proc (x) local i; options operator, arrow; add(p^i*u[i](x), i = 0 .. N) end proc

 

proc (x) local i; options operator, arrow; add(p^i*w[i](x), i = 0 .. N) end proc

(2)

  HPMEq := collect
           ( (1-p)*diff(F(x), x$2)+p*(diff(F(x), x$2)-k1*diff(F(x), x)-k2*F(x)),
              p
           );
  HPMEr := collect
           ( (1-p)*diff(Theta(x),x$2)+p*(diff(Theta(x),x$2)-k11*diff(Theta(x),x)+k12*diff(U(x),x)^2+k13*diff(W(x),x)^2+k14*Theta(x)),
              p
           );
  HPMEs := collect
           ( (1-p)*diff(U(x),x$2)+p*(diff(U(x),x$2)-R*diff(U(x),x)-A-k8*W(x)-k7*U(x)+k5*Theta(x)+k6*F(x)),
              p
           );
  HPMEt := collect
           ( (1-p)*diff(W(x),x$2)+p*(diff(W(x),x$2)-R*diff(W(x),x)+k9*U(x)-k10*W(x)),
              p
           );

(-k1*(diff(f[4](x), x))-k2*f[4](x))*p^5+(-k1*(diff(f[3](x), x))-k2*f[3](x)+diff(diff(f[4](x), x), x))*p^4+(-k1*(diff(f[2](x), x))-k2*f[2](x)+diff(diff(f[3](x), x), x))*p^3+(-k1*(diff(f[1](x), x))-k2*f[1](x)+diff(diff(f[2](x), x), x))*p^2+(-k1*(diff(f[0](x), x))-k2*f[0](x)+diff(diff(f[1](x), x), x))*p+diff(diff(f[0](x), x), x)

 

(k12*(diff(u[4](x), x))^2+k13*(diff(w[4](x), x))^2)*p^9+(2*k12*(diff(u[3](x), x))*(diff(u[4](x), x))+2*k13*(diff(w[3](x), x))*(diff(w[4](x), x)))*p^8+(k12*(2*(diff(u[2](x), x))*(diff(u[4](x), x))+(diff(u[3](x), x))^2)+k13*(2*(diff(w[2](x), x))*(diff(w[4](x), x))+(diff(w[3](x), x))^2))*p^7+(k12*(2*(diff(u[1](x), x))*(diff(u[4](x), x))+2*(diff(u[2](x), x))*(diff(u[3](x), x)))+k13*(2*(diff(w[1](x), x))*(diff(w[4](x), x))+2*(diff(w[2](x), x))*(diff(w[3](x), x))))*p^6+(-k11*(diff(theta[4](x), x))+k12*(2*(diff(u[0](x), x))*(diff(u[4](x), x))+2*(diff(u[1](x), x))*(diff(u[3](x), x))+(diff(u[2](x), x))^2)+k13*(2*(diff(w[0](x), x))*(diff(w[4](x), x))+2*(diff(w[1](x), x))*(diff(w[3](x), x))+(diff(w[2](x), x))^2)+k14*theta[4](x))*p^5+(diff(diff(theta[4](x), x), x)-k11*(diff(theta[3](x), x))+k12*(2*(diff(u[0](x), x))*(diff(u[3](x), x))+2*(diff(u[1](x), x))*(diff(u[2](x), x)))+k13*(2*(diff(w[0](x), x))*(diff(w[3](x), x))+2*(diff(w[1](x), x))*(diff(w[2](x), x)))+k14*theta[3](x))*p^4+(diff(diff(theta[3](x), x), x)-k11*(diff(theta[2](x), x))+k12*(2*(diff(u[0](x), x))*(diff(u[2](x), x))+(diff(u[1](x), x))^2)+k13*(2*(diff(w[0](x), x))*(diff(w[2](x), x))+(diff(w[1](x), x))^2)+k14*theta[2](x))*p^3+(2*k12*(diff(u[0](x), x))*(diff(u[1](x), x))+2*k13*(diff(w[0](x), x))*(diff(w[1](x), x))-k11*(diff(theta[1](x), x))+k14*theta[1](x)+diff(diff(theta[2](x), x), x))*p^2+(k12*(diff(u[0](x), x))^2+k13*(diff(w[0](x), x))^2-k11*(diff(theta[0](x), x))+k14*theta[0](x)+diff(diff(theta[1](x), x), x))*p+diff(diff(theta[0](x), x), x)

 

(-R*(diff(u[4](x), x))-k8*w[4](x)-k7*u[4](x)+k5*theta[4](x)+k6*f[4](x))*p^5+(-R*(diff(u[3](x), x))-k8*w[3](x)-k7*u[3](x)+k5*theta[3](x)+k6*f[3](x)+diff(diff(u[4](x), x), x))*p^4+(-R*(diff(u[2](x), x))-k8*w[2](x)-k7*u[2](x)+k5*theta[2](x)+k6*f[2](x)+diff(diff(u[3](x), x), x))*p^3+(-R*(diff(u[1](x), x))-k8*w[1](x)-k7*u[1](x)+k5*theta[1](x)+k6*f[1](x)+diff(diff(u[2](x), x), x))*p^2+(-R*(diff(u[0](x), x))-k8*w[0](x)-k7*u[0](x)+k5*theta[0](x)+k6*f[0](x)-A+diff(diff(u[1](x), x), x))*p+diff(diff(u[0](x), x), x)

 

(-k10*w[4](x)+k9*u[4](x)-R*(diff(w[4](x), x)))*p^5+(-k10*w[3](x)+k9*u[3](x)-R*(diff(w[3](x), x))+diff(diff(w[4](x), x), x))*p^4+(-k10*w[2](x)+k9*u[2](x)-R*(diff(w[2](x), x))+diff(diff(w[3](x), x), x))*p^3+(-k10*w[1](x)+k9*u[1](x)-R*(diff(w[1](x), x))+diff(diff(w[2](x), x), x))*p^2+(-k10*w[0](x)+k9*u[0](x)-R*(diff(w[0](x), x))+diff(diff(w[1](x), x), x))*p+diff(diff(w[0](x), x), x)

(3)

#
# renamed equ[1] to equ1
#
  for i from 0 to N do
      equ1[i] := coeff(HPMEq, p, i) = 0;
  end do;
#
# renamed equa[1] to equ2
#
  for i from 0 to N do
      equ2[i] := coeff(HPMEr, p, i) = 0;
  end do;
#
# renamed equat[1] to equ3
#
  for i from 0 to N do
      equ3[i] := coeff(HPMEs, p, i) = 0;
  end do;
#
# renamed equati[1] to equ4
#
  for i from 0 to N do
      equ4[i] := coeff(HPMEt, p, i) = 0;
  end do

diff(diff(f[0](x), x), x) = 0

 

-k1*(diff(f[0](x), x))-k2*f[0](x)+diff(diff(f[1](x), x), x) = 0

 

-k1*(diff(f[1](x), x))-k2*f[1](x)+diff(diff(f[2](x), x), x) = 0

 

-k1*(diff(f[2](x), x))-k2*f[2](x)+diff(diff(f[3](x), x), x) = 0

 

-k1*(diff(f[3](x), x))-k2*f[3](x)+diff(diff(f[4](x), x), x) = 0

 

diff(diff(theta[0](x), x), x) = 0

 

k12*(diff(u[0](x), x))^2+k13*(diff(w[0](x), x))^2-k11*(diff(theta[0](x), x))+k14*theta[0](x)+diff(diff(theta[1](x), x), x) = 0

 

2*k12*(diff(u[0](x), x))*(diff(u[1](x), x))+2*k13*(diff(w[0](x), x))*(diff(w[1](x), x))-k11*(diff(theta[1](x), x))+k14*theta[1](x)+diff(diff(theta[2](x), x), x) = 0

 

diff(diff(theta[3](x), x), x)-k11*(diff(theta[2](x), x))+k12*(2*(diff(u[0](x), x))*(diff(u[2](x), x))+(diff(u[1](x), x))^2)+k13*(2*(diff(w[0](x), x))*(diff(w[2](x), x))+(diff(w[1](x), x))^2)+k14*theta[2](x) = 0

 

diff(diff(theta[4](x), x), x)-k11*(diff(theta[3](x), x))+k12*(2*(diff(u[0](x), x))*(diff(u[3](x), x))+2*(diff(u[1](x), x))*(diff(u[2](x), x)))+k13*(2*(diff(w[0](x), x))*(diff(w[3](x), x))+2*(diff(w[1](x), x))*(diff(w[2](x), x)))+k14*theta[3](x) = 0

 

diff(diff(u[0](x), x), x) = 0

 

-R*(diff(u[0](x), x))-k8*w[0](x)-k7*u[0](x)+k5*theta[0](x)+k6*f[0](x)-A+diff(diff(u[1](x), x), x) = 0

 

-R*(diff(u[1](x), x))-k8*w[1](x)-k7*u[1](x)+k5*theta[1](x)+k6*f[1](x)+diff(diff(u[2](x), x), x) = 0

 

-R*(diff(u[2](x), x))-k8*w[2](x)-k7*u[2](x)+k5*theta[2](x)+k6*f[2](x)+diff(diff(u[3](x), x), x) = 0

 

-R*(diff(u[3](x), x))-k8*w[3](x)-k7*u[3](x)+k5*theta[3](x)+k6*f[3](x)+diff(diff(u[4](x), x), x) = 0

 

diff(diff(w[0](x), x), x) = 0

 

-k10*w[0](x)+k9*u[0](x)-R*(diff(w[0](x), x))+diff(diff(w[1](x), x), x) = 0

 

-k10*w[1](x)+k9*u[1](x)-R*(diff(w[1](x), x))+diff(diff(w[2](x), x), x) = 0

 

-k10*w[2](x)+k9*u[2](x)-R*(diff(w[2](x), x))+diff(diff(w[3](x), x), x) = 0

 

-k10*w[3](x)+k9*u[3](x)-R*(diff(w[3](x), x))+diff(diff(w[4](x), x), x) = 0

(4)

#
# Sequentially solve the system(s)
#
  for i from 0 to N do
      cons:= f[i](-1) = 1, f[i](1) = 1,
             theta[i](-1)=1, theta[i](1)=1,
             u[i](-1)=1, u[i](1)=1,
             w[i](-1)=1, w[i](1)=1;
      assign(dsolve( [ equ1[i], equ2[i], equ3[i], equ4[i], cons ]));
  od:

#
# What happens next??
#
# OP can access any of the "individual" functions,
# for example theta[2](x) or w[3](x), just by entering
# these expressions, as in
#
  collect(theta[2](x), x);
  collect(w[3](x), x);
#
# Or one can get the "sum" terms for these functions
# just by entering F(x), Theta(x), U(x) or W(x). These
# expresssions are a bit lengthy
#
# For example, U(x) is shown below
#
  collect(U(x), x);

(1/24)*k14^2*x^4-(1/6)*k14*k11*x^3+(-(1/4)*k14^2-(1/2)*k14)*x^2+(1/6)*k11*k14*x+(5/24)*k14^2+(1/2)*k14+1

 

1+(-(17/180)*A*R*k9+(17/180)*k10^2*R-(17/180)*k10*R*k9+(17/180)*k5*R*k9+(17/180)*R*k9*k6-(17/180)*k7*R*k9-(17/180)*k8*R*k9-(1/6)*k10*R+(1/6)*R*k9)*x+((1/720)*k8*k9^2-(1/720)*k14*k5*k9+(-(1/720)*k7*k8-(1/720)*k7^2+((1/720)*k6+(1/720)*k5-(1/720)*A)*k7+(1/720)*k2*k6)*k9+(-(1/360)*k8-(1/720)*k7+(1/720)*k6+(1/720)*k5-(1/720)*A)*k9*k10-(1/720)*k10^2*k9+(1/720)*k10^3)*x^6+((-(1/60)*k8*R-(1/60)*k7*R+(1/60)*R*k6+(1/60)*k5*R-(1/60)*A*R)*k9-(1/60)*k10*R*k9+(1/60)*k10^2*R)*x^5+(-(1/48)*k8*k9^2+(1/48)*k14*k5*k9+(((1/48)*k7-1/24)*k8+(1/48)*k7^2+(-(1/48)*k6-(1/48)*k5-1/24+(1/48)*A)*k7+(-(1/48)*k2+1/24)*k6+(1/24)*k5-(1/24)*R^2)*k9+(((1/24)*k8+(1/48)*k7-(1/48)*k6-(1/48)*k5-1/24+(1/48)*A)*k9+(1/24)*R^2)*k10+((1/48)*k9+1/24)*k10^2-(1/48)*k10^3)*x^4+(((1/9)*k8*R+(1/9)*k7*R-(1/9)*R*k6-(1/9)*k5*R+((1/9)*A-1/6)*R)*k9+((1/9)*R*k9+(1/6)*R)*k10-(1/9)*k10^2*R)*x^3+((5/48)*k8*k9^2-(5/48)*k14*k5*k9+((-(5/48)*k7+1/4)*k8-(5/48)*k7^2+((5/48)*k6+(5/48)*k5+1/4-(5/48)*A)*k7+((5/48)*k2-1/4)*k6-(1/4)*k5-1/2+(1/12)*R^2)*k9+((-(5/24)*k8-(5/48)*k7+(5/48)*k6+(5/48)*k5+1/4-(5/48)*A)*k9+1/2-(1/12)*R^2)*k10+(-(5/48)*k9-1/4)*k10^2+(5/48)*k10^3)*x^2-(5/24)*k10*k9+(5/24)*k5*k9+(5/24)*k6*k9-(5/24)*k7*k9-(5/24)*k8*k9-(61/720)*k10^3+(5/24)*k10^2+(1/24)*R^2*k10-(1/24)*R^2*k9+(61/720)*k10^2*k9+(61/720)*k7^2*k9-(61/720)*k8*k9^2+(61/720)*A*k10*k9+(61/720)*A*k7*k9-(61/720)*k10*k5*k9-(61/720)*k10*k6*k9+(61/720)*k10*k7*k9+(61/360)*k10*k8*k9+(61/720)*k14*k5*k9-(61/720)*k2*k6*k9-(61/720)*k5*k7*k9-(61/720)*k6*k7*k9+(61/720)*k7*k8*k9+(1/2)*k9-(1/2)*k10

 

p^4*((1/40320)*k8^2*k9^2+(1/40320)*k10^3*k8+(1/40320)*k14^3*k5-(1/40320)*k2^3*k6-(1/40320)*k2^2*k6*k7-(1/40320)*k2*k6*k7^2+(-(1/40320)*k6-(1/40320)*k5+(1/40320)*A)*k7^3+(1/40320)*k7^4+(1/40320)*k7^3*k8+(-(1/20160)*k7*k8^2+(-(1/13440)*k7^2+((1/20160)*k6+(1/20160)*k5-(1/20160)*A)*k7+(1/40320)*k2*k6)*k8)*k9+((-(1/20160)*k8^2+(-(1/20160)*k7+(1/40320)*k6+(1/40320)*k5-(1/40320)*A)*k8)*k9+(1/40320)*k7^2*k8)*k10+(-(1/40320)*k8*k9+(1/40320)*k7*k8)*k10^2+(-(1/40320)*k5*k8*k9+(1/40320)*k5*k7^2)*k14-(1/40320)*k14^2*k5*k7)*x^8+p^4*((-(1/2520)*k1-(1/5040)*R)*k2^2*k6+(-(1/5040)*k1-(1/2520)*R)*k2*k6*k7+(-(1/1680)*R*k6-(1/1680)*k5*R+(1/1680)*A*R)*k7^2+(1/1680)*R*k7^3+(1/1680)*R*k7^2*k8+(-(1/1680)*k8^2*R+(-(1/840)*k7*R+(1/1680)*R*k6+(1/1680)*k5*R-(1/1680)*A*R)*k8)*k9+(-(1/1680)*k8*R*k9+(1/1680)*k7*k8*R)*k10+(1/1680)*R*k10^2*k8+((1/5040)*k11*k5*k7+(1/2520)*k5*k7*R)*k14+(-(1/2520)*k5*k11-(1/5040)*k5*R)*k14^2)*x^7+(p^3*((-(1/720)*k8^2+(-(1/360)*k7+(1/720)*k6+(1/720)*k5-(1/720)*A)*k8)*k9+(1/720)*k7^2*k8+(1/720)*k7^3+(-(1/720)*k6-(1/720)*k5+(1/720)*A)*k7^2-(1/720)*k2*k6*k7-(1/720)*k2^2*k6+(1/720)*k14*k5*k7+(-(1/720)*k8*k9+(1/720)*k7*k8)*k10-(1/720)*k14^2*k5+(1/720)*k10^2*k8)+p^4*(-(1/1440)*k8^2*k9^2-(1/1440)*k10^3*k8-(1/1440)*k14^3*k5+((1/1440)*k2^3-(1/720)*k2^2+(-(1/720)*k1^2-(1/720)*k1*R-(1/720)*R^2)*k2)*k6+(((1/1440)*k2^2-(1/720)*k2-(1/240)*R^2)*k6-(1/240)*R^2*k5+(1/240)*A*R^2)*k7+(((1/1440)*k2-1/720)*k6-(1/720)*k5+(1/240)*R^2)*k7^2+((1/1440)*k6+(1/1440)*k5-(1/1440)*A+1/720)*k7^3-(1/1440)*k7^4+(-(1/1440)*k7^3+(1/720)*k7^2+(1/240)*R^2*k7)*k8+(((1/720)*k7-1/720)*k8^2+((1/480)*k7^2+(-(1/720)*k6-(1/720)*k5+(1/720)*A-1/360)*k7+(-(1/1440)*k2+1/720)*k6+(1/720)*k5-(1/240)*R^2)*k8)*k9+(((1/720)*k8^2+((1/720)*k7-(1/1440)*k6-(1/1440)*k5+(1/1440)*A-1/720)*k8)*k9+(-(1/1440)*k7^2+(1/720)*k7+(1/240)*R^2)*k8)*k10+((1/1440)*k8*k9+(-(1/1440)*k7+1/720)*k8)*k10^2+((1/360)*k5*k8^2+((1/180)*k5*k7-(1/180)*k5*k6-(1/180)*k5^2+(1/180)*A*k5)*k8+(1/360)*k5*k7^2+(-(1/180)*k5*k6-(1/180)*k5^2+(1/180)*A*k5)*k7+(1/360)*k5*k6^2+((1/180)*k5^2-(1/180)*A*k5)*k6+(1/360)*k5^3-(1/180)*A*k5^2+(1/360)*A^2*k5)*k12+((1/360)*k5*k10^2-(1/180)*k10*k5*k9+(1/360)*k5*k9^2)*k13+((1/720)*k5*k11^2+(1/720)*k11*k5*R+(1/1440)*k5*k8*k9-(1/1440)*k5*k7^2+(1/720)*k5*k7+(1/720)*R^2*k5)*k14+((1/1440)*k5*k7-(1/720)*k5)*k14^2))*x^6+(p^3*(-(1/60)*k8*R*k9+(1/60)*k7*k8*R+(1/60)*k7^2*R+(-(1/60)*R*k6-(1/60)*k5*R+(1/60)*A*R)*k7+(-(1/120)*k1-(1/120)*R)*k2*k6+((1/120)*k5*k11+(1/120)*k5*R)*k14+(1/60)*k10*k8*R)+p^4*((((1/180)*k1+(1/240)*R)*k2^2+(-(1/120)*k1-(1/120)*R)*k2-(1/120)*R^3)*k6-(1/120)*R^3*k5+((((1/720)*k1+(1/120)*R)*k2-(1/60)*R)*k6-(1/60)*k5*R+(1/120)*R^3)*k7+((7/720)*R*k6+(7/720)*k5*R+(1/60-(7/720)*A)*R)*k7^2-(7/720)*R*k7^3+(-(7/720)*k7^2*R+(1/60)*k7*R+(1/120)*R^3)*k8+((7/720)*k8^2*R+((7/360)*k7*R-(7/720)*R*k6-(7/720)*k5*R+(-1/60+(7/720)*A)*R)*k8)*k9+((7/720)*k8*R*k9+(-(7/720)*k7*R+(1/60)*R)*k8)*k10-(7/720)*R*k10^2*k8+((-(1/720)*k5*k7+(1/120)*k5)*k11-(1/120)*k5*k7*R+(1/120)*k5*R)*k14+((1/180)*k5*k11+(1/240)*k5*R)*k14^2+(1/120)*A*R^3))*x^5+(p^2*((1/24)*A*k7+(1/24)*k10*k8+(1/24)*k14*k5-(1/24)*k2*k6-(1/24)*k5*k7-(1/24)*k6*k7+(1/24)*k7^2+(1/24)*k7*k8-(1/24)*k8*k9)+p^3*(((1/48)*k8^2+((1/24)*k7-(1/48)*k6-(1/48)*k5-1/24+(1/48)*A)*k8)*k9+(-(1/48)*k7^2+(1/24)*k7+(1/24)*R^2)*k8-(1/48)*k7^3+((1/48)*k6+(1/48)*k5+1/24-(1/48)*A)*k7^2+(((1/48)*k2-1/24)*k6-(1/24)*k5+(1/24)*R^2)*k7+((1/48)*k2^2-(1/24)*k2-(1/24)*R^2)*k6-(1/24)*R^2*k5+(1/24)*A*R^2+(-(1/48)*k5*k7+(1/24)*k5)*k14+((1/48)*k8*k9+(-(1/48)*k7+1/24)*k8)*k10+(1/48)*k14^2*k5-(1/48)*k10^2*k8)+p^4*((5/576)*k8^2*k9^2+(5/576)*k10^3*k8+(5/576)*k14^3*k5+(-(5/576)*k2^3+(1/48)*k2^2+((1/144)*k1^2+(1/144)*k1*R+(1/48)*R^2-1/24)*k2-(1/24)*R^2)*k6-(1/24)*R^2*k5+((-(5/576)*k2^2+(1/48)*k2+(5/144)*R^2-1/24)*k6+((5/144)*R^2-1/24)*k5+(1/24-(5/144)*A)*R^2)*k7+((-(5/576)*k2+1/48)*k6+(1/48)*k5+1/24-(5/144)*R^2)*k7^2+(-(5/576)*k6-(5/576)*k5+(5/576)*A-1/48)*k7^3+(5/576)*k7^4+((5/576)*k7^3-(1/48)*k7^2+(1/24-(5/144)*R^2)*k7+(1/24)*R^2)*k8+((-(5/288)*k7+1/48)*k8^2+(-(5/192)*k7^2+((5/288)*k6+(5/288)*k5-(5/288)*A+1/24)*k7+((5/576)*k2-1/48)*k6-(1/48)*k5+(5/144)*R^2-1/24)*k8)*k9+((-(5/288)*k8^2+(-(5/288)*k7+(5/576)*k6+(5/576)*k5-(5/576)*A+1/48)*k8)*k9+((5/576)*k7^2-(1/48)*k7+1/24-(5/144)*R^2)*k8)*k10+(-(5/576)*k8*k9+((5/576)*k7-1/48)*k8)*k10^2+(-(1/144)*k5*k11^2-(1/144)*k11*k5*R-(5/576)*k5*k8*k9+(5/576)*k5*k7^2-(1/48)*k5*k7+(1/24-(1/48)*R^2)*k5)*k14+(-(5/576)*k5*k7+(1/48)*k5)*k14^2))*x^4+(p^2*((1/6)*A*R-(1/6)*k5*R-(1/6)*R*k6+(1/6)*k7*R+(1/6)*k8*R)+p^3*((1/9)*k8*R*k9+(-(1/9)*k7*R+(1/6)*R)*k8-(1/9)*k7^2*R+((1/9)*R*k6+(1/9)*k5*R+(1/6-(1/9)*A)*R)*k7+(((1/36)*k1+(1/12)*R)*k2-(1/6)*R)*k6-(1/6)*k5*R+(-(1/36)*k5*k11-(1/12)*k5*R)*k14-(1/9)*k10*k8*R)+p^4*(((-(17/1080)*k1-(5/144)*R)*k2^2+((1/36)*k1+(1/12)*R)*k2+(1/36)*R^3-(1/6)*R)*k6+((1/36)*R^3-(1/6)*R)*k5+(((-(7/2160)*k1-(17/360)*R)*k2+(1/9)*R)*k6+(1/9)*k5*R-(1/36)*R^3+(1/6)*R)*k7+(-(109/2160)*R*k6-(109/2160)*k5*R+((109/2160)*A-1/9)*R)*k7^2+(109/2160)*R*k7^3+((109/2160)*k7^2*R-(1/9)*k7*R-(1/36)*R^3+(1/6)*R)*k8+(-(109/2160)*k8^2*R+(-(109/1080)*k7*R+(109/2160)*R*k6+(109/2160)*k5*R+(-(109/2160)*A+1/9)*R)*k8)*k9+(-(109/2160)*k8*R*k9+((109/2160)*k7*R-(1/9)*R)*k8)*k10+(109/2160)*R*k10^2*k8+(((7/2160)*k5*k7-(1/36)*k5)*k11+(17/360)*k5*k7*R-(1/12)*k5*R)*k14+(-(17/1080)*k5*k11-(5/144)*k5*R)*k14^2-(1/36)*A*R^3))*x^3+(p*((1/2)*k8+(1/2)*k7-(1/2)*k5-(1/2)*k6+(1/2)*A)+p^2*(-(1/4)*A*k7-(1/4)*k10*k8-(1/4)*k14*k5+(1/4)*k2*k6+(1/4)*k5*k7+(1/4)*k6*k7-(1/4)*k7^2-(1/4)*k7*k8+(1/4)*k8*k9-(1/2)*k5-(1/2)*k6+(1/2)*k7+(1/2)*k8)+p^3*((-(5/48)*k8^2+(-(5/24)*k7+(5/48)*k6+(5/48)*k5+1/4-(5/48)*A)*k8)*k9+((5/48)*k7^2-(1/4)*k7+1/2-(1/12)*R^2)*k8+(5/48)*k7^3+(-(5/48)*k6-(5/48)*k5-1/4+(5/48)*A)*k7^2+((-(5/48)*k2+1/4)*k6+(1/4)*k5+1/2-(1/12)*R^2)*k7+(-(5/48)*k2^2+(1/4)*k2-1/2+(1/12)*R^2)*k6+(-1/2+(1/12)*R^2)*k5-(1/12)*A*R^2+((5/48)*k5*k7-(1/4)*k5)*k14+(-(5/48)*k8*k9+((5/48)*k7-1/4)*k8)*k10-(5/48)*k14^2*k5+(5/48)*k10^2*k8)+p^4*(-(61/1440)*k8^2*k9^2-(61/1440)*k10^3*k8-(61/1440)*k14^3*k5+((61/1440)*k2^3-(5/48)*k2^2+(-(1/48)*k1^2-(7/720)*k1*R-(3/80)*R^2+1/4)*k2-1/2+(1/12)*R^2)*k6+(-1/2+(1/12)*R^2)*k5+(((61/1440)*k2^2-(5/48)*k2+1/4-(49/720)*R^2)*k6+(1/4-(49/720)*R^2)*k5+1/2+((49/720)*A-1/12)*R^2)*k7+(((61/1440)*k2-5/48)*k6-(5/48)*k5-1/4+(49/720)*R^2)*k7^2+((61/1440)*k6+(61/1440)*k5-(61/1440)*A+5/48)*k7^3-(61/1440)*k7^4+(-(61/1440)*k7^3+(5/48)*k7^2+(-1/4+(49/720)*R^2)*k7+1/2-(1/12)*R^2)*k8+(((61/720)*k7-5/48)*k8^2+((61/480)*k7^2+(-(61/720)*k6-(61/720)*k5-5/24+(61/720)*A)*k7+(-(61/1440)*k2+5/48)*k6+(5/48)*k5+1/4-(49/720)*R^2)*k8)*k9+(((61/720)*k8^2+((61/720)*k7-(61/1440)*k6-(61/1440)*k5+(61/1440)*A-5/48)*k8)*k9+(-(61/1440)*k7^2+(5/48)*k7-1/4+(49/720)*R^2)*k8)*k10+((61/1440)*k8*k9+(-(61/1440)*k7+5/48)*k8)*k10^2+(-(1/24)*k5*k8^2+(-(1/12)*k5*k7+(1/12)*k5*k6+(1/12)*k5^2-(1/12)*A*k5)*k8-(1/24)*k5*k7^2+((1/12)*k5*k6+(1/12)*k5^2-(1/12)*A*k5)*k7-(1/24)*k5*k6^2+(-(1/12)*k5^2+(1/12)*A*k5)*k6-(1/24)*k5^3+(1/12)*A*k5^2-(1/24)*A^2*k5)*k12+(-(1/24)*k5*k10^2+(1/12)*k10*k5*k9-(1/24)*k5*k9^2)*k13+((1/48)*k5*k11^2+(7/720)*k11*k5*R+(61/1440)*k5*k8*k9-(61/1440)*k5*k7^2+(5/48)*k5*k7+(-1/4+(3/80)*R^2)*k5)*k14+((61/1440)*k5*k7-(5/48)*k5)*k14^2))*x^2+(p^2*(-(1/6)*A*R+(1/6)*k5*R+(1/6)*R*k6-(1/6)*k7*R-(1/6)*k8*R)+p^3*((17/180)*A*k7*R+(17/180)*k10*k8*R+(3/40)*k5*R*k14-(3/40)*k2*R*k6-(17/180)*k5*k7*R-(17/180)*k7*R*k6+(17/180)*k7^2*R+(17/180)*k7*k8*R-(17/180)*k8*R*k9-(7/360)*k2*k1*k6+(7/360)*k5*k11*k14+(1/6)*k5*R+(1/6)*R*k6-(1/6)*k7*R-(1/6)*k8*R)+p^4*((125/3024)*A*R*k8*k9-(125/3024)*R*k10*k7*k8-(11/280)*R*k14*k5*k7+(11/280)*R*k2*k6*k7-(125/3024)*R*k5*k8*k9-(125/3024)*R*k6*k8*k9+(31/15120)*k1*k2*k6*k7-(31/15120)*k11*k14*k5*k7-(1/6)*k8*R-(1/6)*k7*R+(1/6)*k5*R+(1/6)*R*k6-(17/180)*k8*R*k9+(125/3024)*R*k10*k8*k9+(125/1512)*R*k7*k8*k9+(17/180)*k7^2*R+(17/180)*k10*k8*R+(3/40)*k5*R*k14-(3/40)*k2*R*k6-(17/180)*k5*k7*R-(17/180)*k7*R*k6+(17/180)*k7*k8*R-(7/360)*k2*k1*k6+(7/360)*k5*k11*k14+(7/360)*A*R^3-(7/360)*R^3*k5-(7/360)*R^3*k6+(7/360)*R^3*k7+(7/360)*R^3*k8-(125/3024)*R*k7^3-(125/3024)*A*R*k7^2-(125/3024)*R*k10^2*k8+(31/1008)*R*k14^2*k5+(31/1008)*R*k2^2*k6+(125/3024)*R*k5*k7^2+(125/3024)*R*k6*k7^2-(125/3024)*R*k7^2*k8+(125/3024)*R*k8^2*k9+(2/189)*k1*k2^2*k6+(2/189)*k11*k14^2*k5))*x+1+p*(-(1/2)*k8-(1/2)*k7+(1/2)*k5+(1/2)*k6-(1/2)*A+1)+p^2*((5/24)*A*k7+(5/24)*k10*k8+(5/24)*k14*k5-(5/24)*k2*k6-(5/24)*k5*k7-(5/24)*k6*k7+(5/24)*k7^2+(5/24)*k7*k8-(5/24)*k8*k9+(1/2)*k5+(1/2)*k6-(1/2)*k7-(1/2)*k8+1)+p^3*(1-(1/2)*k8-(1/2)*k7+(1/2)*k5+(1/2)*k6-(5/24)*k8*k9+(5/24)*k10*k8+(5/24)*k14*k5-(5/24)*k2*k6-(5/24)*k5*k7-(5/24)*k6*k7+(5/24)*k7*k8-(61/720)*k7^3+(5/24)*k7^2+(1/24)*A*R^2-(61/720)*A*k7^2-(1/24)*R^2*k5-(1/24)*R^2*k6+(1/24)*R^2*k7+(1/24)*R^2*k8-(61/720)*k10^2*k8+(61/720)*k14^2*k5+(61/720)*k2^2*k6+(61/720)*k5*k7^2+(61/720)*k6*k7^2-(61/720)*k7^2*k8+(61/720)*k8^2*k9+(61/720)*k10*k8*k9+(61/360)*k7*k8*k9+(61/720)*A*k8*k9-(61/720)*k10*k7*k8-(61/720)*k14*k5*k7+(61/720)*k2*k6*k7-(61/720)*k5*k8*k9-(61/720)*k6*k8*k9)+p^4*(1-(1/2)*k8-(1/2)*k7+(1/2)*k5+(1/2)*k6-(5/24)*k8*k9+(5/24)*k10*k8+(5/24)*k14*k5-(5/24)*k2*k6-(5/24)*k5*k7-(5/24)*k6*k7+(5/24)*k7*k8+(3/80)*R^2*k8*k9-(277/8064)*k10^2*k8*k9-(277/2688)*k7^2*k8*k9+(7/180)*A^2*k12*k5-(3/80)*A*R^2*k7-(7/90)*A*k12*k5^2-(3/80)*R^2*k10*k8-(13/720)*R^2*k14*k5+(13/720)*R^2*k2*k6+(3/80)*R^2*k5*k7+(3/80)*R^2*k6*k7-(3/80)*R^2*k7*k8+(11/720)*k1^2*k2*k6+(7/180)*k10^2*k13*k5+(277/8064)*k10^2*k7*k8+(277/8064)*k10*k7^2*k8-(277/4032)*k10*k8^2*k9-(11/720)*k11^2*k14*k5+(7/90)*k12*k5^2*k6-(7/90)*k12*k5^2*k7-(7/90)*k12*k5^2*k8+(7/180)*k12*k5*k6^2+(7/180)*k12*k5*k7^2+(7/180)*k12*k5*k8^2+(7/180)*k13*k5*k9^2-(277/8064)*k14^2*k5*k7+(277/8064)*k14*k5*k7^2-(277/8064)*k2^2*k6*k7-(277/8064)*k2*k6*k7^2-(277/4032)*k7*k8^2*k9-(61/720)*k7^3+(5/24)*k7^2+(277/8064)*A*k7^3-(3/80)*R^2*k7^2+(277/8064)*k10^3*k8+(7/180)*k12*k5^3+(277/8064)*k14^3*k5-(277/8064)*k2^3*k6-(277/8064)*k5*k7^3-(277/8064)*k6*k7^3+(277/8064)*k7^3*k8+(277/8064)*k8^2*k9^2-(277/4032)*k10*k7*k8*k9-(1/24)*R^2*k5-(1/24)*R^2*k6+(1/24)*R^2*k7+(1/24)*R^2*k8-(61/720)*k10^2*k8+(61/720)*k14^2*k5+(61/720)*k2^2*k6+(61/720)*k5*k7^2+(61/720)*k6*k7^2-(61/720)*k7^2*k8+(61/720)*k8^2*k9+(61/720)*k10*k8*k9+(61/360)*k7*k8*k9-(61/720)*k10*k7*k8-(61/720)*k14*k5*k7+(61/720)*k2*k6*k7-(61/720)*k5*k8*k9-(61/720)*k6*k8*k9-(277/8064)*A*k10*k8*k9-(7/90)*A*k12*k5*k6+(7/90)*A*k12*k5*k7+(7/90)*A*k12*k5*k8-(277/4032)*A*k7*k8*k9+(1/240)*R*k1*k2*k6-(1/240)*R*k11*k14*k5-(7/90)*k10*k13*k5*k9+(277/8064)*k10*k5*k8*k9+(277/8064)*k10*k6*k8*k9-(7/90)*k12*k5*k6*k7-(7/90)*k12*k5*k6*k8+(7/90)*k12*k5*k7*k8-(277/8064)*k14*k5*k8*k9+(277/8064)*k2*k6*k8*k9+(277/4032)*k5*k7*k8*k9+(277/4032)*k6*k7*k8*k9+(277/8064)*k7^4)

(5)

 


 

Download HOM.mw

 

Defining a function with the '->' notation only really allows very simple functions. In particular 'if' statements are only allowed to be of the form

if(condition, doThis, otherwiseDoThis)

'elif' clauses are not permitted. Thus your definiton of the function 'q1' is invalid. Interestingly, you have perfectly valid definition of 'q1' as a proc(), which is commented out! If I comment out the 'q1->....' statement and uncomment the qi:=proc(....) statement, then all I get are a few warnings in the calling procedure about undeclared local variables. This is easily fixed by adding these to the 'local' statement.

The result executes with no errors - obviously I have not checked whether the resulting procedure provide the answers you want, bu the attached executes with no errors.

restart; with(Physics); with(ExcelTools); with(LinearAlgebra); with(Statistics); with(PolynomialTools); with(ArrayTools); with(ListTools); with(DocumentTools); with(combinat); with(linalg); with(plots); interface(rtablesize = 50)

Coupled_Freq_8th := proc (freq, terms, kmax, jmax, imax) local w1, w2, w3, q0, q1, q2, q3, q4, q5, q6, q7, q8, q9, q10, q11, q12, n1, n2, n3, nn, largest, nq, nqn, q1mat, q2mat, q3mat, q4mat, q5mat, q6mat, q7mat, q8mat, q9mat, q10mat, allqmat, h_harm, h_q1, h_q2, h_q3, h_q12_m3, h_q12_m4, h_q12_m5, h_q12_m6, h_q12_m7, h_q12_m8, h_q12_m9, h_q12_m10, h_q12, h_q13_m3, h_q13_m4, h_q13_m5, h_q13_m6, h_q13_m7, h_q13_m8, h_q13_m9, h_q13_m10, h_q13, h_q23_m3, h_q23_m4, h_q23_m5, h_q23_m6, h_q23_m7, h_q23_m8, h_q23_m9, h_q23_m10, h_q23, evals_no, v1no, v2no, v3no, uncoupled_freq, evals_pairmat, evals_pairlist, temp, i, j, m, v1, v2, v3, coupled_freq, row1, row2, row3, h_no, h_tot, nq1, nq2, nq12, temp1, eval_size, eval_cont, new4, z, k, v1_t, v2_t, v3_t, coupled_freq_t; w1 := freq[1]; w2 := freq[2]; w3 := freq[3]; q0 := proc (n1, n2) options operator, arrow; `if`(n1 = n2, 1, 0) end proc; q1 := proc (n1, n2) if n1 < 0 then 0 elif n2 < 0 then 0 elif abs(n1+Physics:-`*`(-1, n2)) = 1 then sqrt(Physics:-`*`(max(n1, n2), Physics:-`^`(2, -1))) else 0 end if end proc; q2 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q1(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q3 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q2(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q4 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q3(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q5 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q4(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q6 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q5(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q7 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q6(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q8 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q7(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q9 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q8(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q10 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q9(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q11 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q10(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; q12 := proc (n1, n2) options operator, arrow; add(Physics:-`*`(q11(n1, i), q1(i, n2)), i = n2-1 .. n2+1) end proc; seq(seq(assign(terms[z][i] = terms[z+6][i]), i = 1 .. nops(terms[z])), z = 1 .. 6); n1 := [seq(seq(seq(k, k = 0 .. kmax), j = 0 .. jmax), i = 0 .. imax)]; nn := nops(n1); n2 := [seq(seq(seq(j, k = 0 .. kmax), j = 0 .. jmax), i = 0 .. imax)]; n3 := [seq(seq(seq(i, k = 0 .. kmax), j = 0 .. jmax), i = 0 .. imax)]; row1 := [n1[1], n2[1], n3[1]]; row2 := [n1[2], n2[2], n3[2]]; row3 := [n1[3], n2[3], n3[3]]; largest := max(imax, jmax, kmax); nq := [seq(0 .. largest)]; nqn := nops(nq); q1mat := Matrix(nqn, nqn, [seq(seq(q1(nq[i], nq[j]), j = 1 .. nqn), i = 1 .. nqn)]); q2mat := Matrix(nqn, nqn, [seq(seq(q2(nq[i], nq[j]), j = 1 .. nqn), i = 1 .. nqn)]); q3mat := Matrix(nqn, nqn, [seq(seq(q3(nq[i], nq[j]), j = 1 .. nqn), i = 1 .. nqn)]); q4mat := Matrix(nqn, nqn, [seq(seq(q4(nq[i], nq[j]), j = 1 .. nqn), i = 1 .. nqn)]); q5mat := Matrix(nqn, nqn, [seq(seq(q5(nq[i], nq[j]), j = 1 .. nqn), i = 1 .. nqn)]); q6mat := Matrix(nqn, nqn, [seq(seq(q6(nq[i], nq[j]), j = 1 .. nqn), i = 1 .. nqn)]); q7mat := Matrix(nqn, nqn, [seq(seq(q7(nq[i], nq[j]), j = 1 .. nqn), i = 1 .. nqn)]); q8mat := Matrix(nqn, nqn, [seq(seq(q8(nq[i], nq[j]), j = 1 .. nqn), i = 1 .. nqn)]); allqmat := [q1mat, q2mat, q3mat, q4mat, q5mat, q6mat, q7mat, q8mat]; h_harm := Matrix(nn, nn, [seq(seq(Physics:-`*`(Physics:-`*`(Physics:-`*`(q0(n1[i], n1[j]), q0(n2[i], n2[j])), q0(n3[i], n3[j])), Physics:-`*`(w1, n1[i]+.5)+Physics:-`*`(w2, n2[i]+.5)+Physics:-`*`(w3, n3[i]+.5)), j = 1 .. nn), i = 1 .. nn)]); h_q1 := `~`[Physics:-`*`](Matrix(nn, nn, [seq(seq(Physics:-`*`(Physics:-`*`(q0(n2[i], n2[j]), q0(n3[i], n3[j])), Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], k300)+Physics:-`*`(q4mat[n1[i]+1, n1[j]+1], k400)+Physics:-`*`(q5mat[n1[i]+1, n1[j]+1], k500)+Physics:-`*`(q6mat[n1[i]+1, n1[j]+1], k600)+Physics:-`*`(q7mat[n1[i]+1, n1[j]+1], k700)+Physics:-`*`(q8mat[n1[i]+1, n1[j]+1], k800)), j = 1 .. nn), i = 1 .. nn)]), 219474.6); h_q2 := `~`[Physics:-`*`](Matrix(nn, nn, [seq(seq(Physics:-`*`(Physics:-`*`(q0(n1[i], n1[j]), q0(n3[i], n3[j])), Physics:-`*`(q3mat[n2[i]+1, n2[j]+1], k030)+Physics:-`*`(q4mat[n2[i]+1, n2[j]+1], k040)+Physics:-`*`(q5mat[n2[i]+1, n2[j]+1], k050)+Physics:-`*`(q6mat[n2[i]+1, n2[j]+1], k060)+Physics:-`*`(q7mat[n2[i]+1, n2[j]+1], k070)+Physics:-`*`(q8mat[n2[i]+1, n2[j]+1], k080)), j = 1 .. nn), i = 1 .. nn)]), 219474.6); h_q3 := `~`[Physics:-`*`](Matrix(nn, nn, [seq(seq(Physics:-`*`(Physics:-`*`(q0(n1[i], n1[j]), q0(n2[i], n2[j])), Physics:-`*`(q3mat[n3[i]+1, n3[j]+1], k003)+Physics:-`*`(q4mat[n3[i]+1, n3[j]+1], k004)+Physics:-`*`(q5mat[n3[i]+1, n3[j]+1], k005)+Physics:-`*`(q6mat[n3[i]+1, n3[j]+1], k006)+Physics:-`*`(q7mat[n3[i]+1, n3[j]+1], k007)+Physics:-`*`(q8mat[n3[i]+1, n3[j]+1], k008)), j = 1 .. nn), i = 1 .. nn)]), 219474.6); h_q12_m3 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n3[i], n3[j]), Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q1mat[n2[i]+1, n2[j]+1]), k210)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q2mat[n2[i]+1, n2[j]+1]), k120)), j = 1 .. nn), i = 1 .. nn)]); h_q12_m4 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n3[i], n3[j]), Physics:-`*`(Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], q1mat[n2[i]+1, n2[j]+1]), k310)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q2mat[n2[i]+1, n2[j]+1]), k220)), j = 1 .. nn), i = 1 .. nn)]); h_q12_m5 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n3[i], n3[j]), Physics:-`*`(Physics:-`*`(q4mat[n1[i]+1, n1[j]+1], q1mat[n2[i]+1, n2[j]+1]), k410)+Physics:-`*`(Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], q2mat[n2[i]+1, n2[j]+1]), k320)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q3mat[n2[i]+1, n2[j]+1]), k230)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q4mat[n2[i]+1, n2[j]+1]), k140)), j = 1 .. nn), i = 1 .. nn)]); h_q12_m6 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n3[i], n3[j]), Physics:-`*`(Physics:-`*`(q5mat[n1[i]+1, n1[j]+1], q1mat[n2[i]+1, n2[j]+1]), k510)+Physics:-`*`(Physics:-`*`(q4mat[n1[i]+1, n1[j]+1], q2mat[n2[i]+1, n2[j]+1]), k420)+Physics:-`*`(Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], q3mat[n2[i]+1, n2[j]+1]), k330)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q4mat[n2[i]+1, n2[j]+1]), k240)), j = 1 .. nn), i = 1 .. nn)]); h_q12_m7 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n3[i], n3[j]), Physics:-`*`(Physics:-`*`(q6mat[n1[i]+1, n1[j]+1], q1mat[n2[i]+1, n2[j]+1]), k610)+Physics:-`*`(Physics:-`*`(q5mat[n1[i]+1, n1[j]+1], q2mat[n2[i]+1, n2[j]+1]), k520)+Physics:-`*`(Physics:-`*`(q4mat[n1[i]+1, n1[j]+1], q3mat[n2[i]+1, n2[j]+1]), k430)+Physics:-`*`(Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], q4mat[n2[i]+1, n2[j]+1]), k340)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q5mat[n2[i]+1, n2[j]+1]), k250)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q6mat[n2[i]+1, n2[j]+1]), k160)), j = 1 .. nn), i = 1 .. nn)]); h_q12_m8 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n3[i], n3[j]), Physics:-`*`(Physics:-`*`(q7mat[n1[i]+1, n1[j]+1], q1mat[n2[i]+1, n2[j]+1]), k710)+Physics:-`*`(Physics:-`*`(q6mat[n1[i]+1, n1[j]+1], q2mat[n2[i]+1, n2[j]+1]), k620)+Physics:-`*`(Physics:-`*`(q5mat[n1[i]+1, n1[j]+1], q3mat[n2[i]+1, n2[j]+1]), k530)+Physics:-`*`(Physics:-`*`(q4mat[n1[i]+1, n1[j]+1], q4mat[n2[i]+1, n2[j]+1]), k440)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q6mat[n2[i]+1, n2[j]+1]), k260)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q7mat[n2[i]+1, n2[j]+1]), k170)), j = 1 .. nn), i = 1 .. nn)]); h_q12 := `~`[Physics:-`*`](h_q12_m3+h_q12_m4+h_q12_m5+h_q12_m6+h_q12_m7+h_q12_m8, 219474.6); h_q13_m3 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n2[i], n2[j]), Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q1mat[n3[i]+1, n3[j]+1]), k201)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q2mat[n3[i]+1, n3[j]+1]), k102)), j = 1 .. nn), i = 1 .. nn)]); h_q13_m4 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n2[i], n2[j]), Physics:-`*`(Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], q1mat[n3[i]+1, n3[j]+1]), k301)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q2mat[n3[i]+1, n3[j]+1]), k202)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q3mat[n3[i]+1, n3[j]+1]), k103)), j = 1 .. nn), i = 1 .. nn)]); h_q13_m5 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n2[i], n2[j]), Physics:-`*`(Physics:-`*`(q4mat[n1[i]+1, n1[j]+1], q1mat[n3[i]+1, n3[j]+1]), k401)+Physics:-`*`(Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], q2mat[n3[i]+1, n3[j]+1]), k302)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q3mat[n3[i]+1, n3[j]+1]), k203)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q4mat[n3[i]+1, n3[j]+1]), k104)), j = 1 .. nn), i = 1 .. nn)]); h_q13_m6 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n2[i], n2[j]), Physics:-`*`(Physics:-`*`(q5mat[n1[i]+1, n1[j]+1], q1mat[n3[i]+1, n3[j]+1]), k501)+Physics:-`*`(Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], q3mat[n3[i]+1, n3[j]+1]), k303)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q4mat[n3[i]+1, n3[j]+1]), k204)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q5mat[n3[i]+1, n3[j]+1]), k105)), j = 1 .. nn), i = 1 .. nn)]); h_q13_m7 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n2[i], n2[j]), Physics:-`*`(Physics:-`*`(q6mat[n1[i]+1, n1[j]+1], q1mat[n3[i]+1, n3[j]+1]), k601)+Physics:-`*`(Physics:-`*`(q5mat[n1[i]+1, n1[j]+1], q2mat[n3[i]+1, n3[j]+1]), k502)+Physics:-`*`(Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], q4mat[n3[i]+1, n3[j]+1]), k304)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q5mat[n3[i]+1, n3[j]+1]), k205)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q6mat[n3[i]+1, n3[j]+1]), k106)), j = 1 .. nn), i = 1 .. nn)]); h_q13_m8 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n2[i], n2[j]), Physics:-`*`(Physics:-`*`(q7mat[n1[i]+1, n1[j]+1], q1mat[n3[i]+1, n3[j]+1]), k701)+Physics:-`*`(Physics:-`*`(q6mat[n1[i]+1, n1[j]+1], q2mat[n3[i]+1, n3[j]+1]), k602)+Physics:-`*`(Physics:-`*`(q5mat[n1[i]+1, n1[j]+1], q3mat[n3[i]+1, n3[j]+1]), k503)+Physics:-`*`(Physics:-`*`(q4mat[n1[i]+1, n1[j]+1], q4mat[n3[i]+1, n3[j]+1]), k404)+Physics:-`*`(Physics:-`*`(q3mat[n1[i]+1, n1[j]+1], q5mat[n3[i]+1, n3[j]+1]), k305)+Physics:-`*`(Physics:-`*`(q2mat[n1[i]+1, n1[j]+1], q6mat[n3[i]+1, n3[j]+1]), k206)+Physics:-`*`(Physics:-`*`(q1mat[n1[i]+1, n1[j]+1], q7mat[n3[i]+1, n3[j]+1]), k107)), j = 1 .. nn), i = 1 .. nn)]); h_q13 := `~`[Physics:-`*`](h_q13_m3+h_q13_m4+h_q13_m5+h_q13_m6+h_q13_m7+h_q13_m8, 219474.6); h_q23_m3 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n1[i], n1[j]), Physics:-`*`(Physics:-`*`(q2mat[n2[i]+1, n2[j]+1], q1mat[n3[i]+1, n3[j]+1]), k021)+Physics:-`*`(Physics:-`*`(q1mat[n2[i]+1, n2[j]+1], q2mat[n3[i]+1, n3[j]+1]), k012)), j = 1 .. nn), i = 1 .. nn)]); h_q23_m4 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n1[i], n1[j]), Physics:-`*`(Physics:-`*`(q3mat[n2[i]+1, n2[j]+1], q1mat[n3[i]+1, n3[j]+1]), k031)+Physics:-`*`(Physics:-`*`(q2mat[n2[i]+1, n2[j]+1], q2mat[n3[i]+1, n3[j]+1]), k022)+Physics:-`*`(Physics:-`*`(q1mat[n2[i]+1, n2[j]+1], q3mat[n3[i]+1, n3[j]+1]), k013)), j = 1 .. nn), i = 1 .. nn)]); h_q23_m5 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n1[i], n1[j]), Physics:-`*`(Physics:-`*`(q4mat[n2[i]+1, n2[j]+1], q1mat[n3[i]+1, n3[j]+1]), k041)+Physics:-`*`(Physics:-`*`(q3mat[n2[i]+1, n2[j]+1], q2mat[n3[i]+1, n3[j]+1]), k032)+Physics:-`*`(Physics:-`*`(q2mat[n2[i]+1, n2[j]+1], q3mat[n3[i]+1, n3[j]+1]), k023)+Physics:-`*`(Physics:-`*`(q1mat[n2[i]+1, n2[j]+1], q4mat[n3[i]+1, n3[j]+1]), k014)), j = 1 .. nn), i = 1 .. nn)]); h_q23_m6 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n1[i], n1[j]), Physics:-`*`(Physics:-`*`(q5mat[n2[i]+1, n2[j]+1], q1mat[n3[i]+1, n3[j]+1]), k051)+Physics:-`*`(Physics:-`*`(q4mat[n2[i]+1, n2[j]+1], q2mat[n3[i]+1, n3[j]+1]), k042)+Physics:-`*`(Physics:-`*`(q3mat[n2[i]+1, n2[j]+1], q3mat[n3[i]+1, n3[j]+1]), k033)+Physics:-`*`(Physics:-`*`(q2mat[n2[i]+1, n2[j]+1], q4mat[n3[i]+1, n3[j]+1]), k024)+Physics:-`*`(Physics:-`*`(q1mat[n2[i]+1, n2[j]+1], q5mat[n3[i]+1, n3[j]+1]), k015)), j = 1 .. nn), i = 1 .. nn)]); h_q23_m7 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n1[i], n1[j]), Physics:-`*`(Physics:-`*`(q6mat[n2[i]+1, n2[j]+1], q1mat[n3[i]+1, n3[j]+1]), k061)+Physics:-`*`(Physics:-`*`(q5mat[n2[i]+1, n2[j]+1], q2mat[n3[i]+1, n3[j]+1]), k052)+Physics:-`*`(Physics:-`*`(q4mat[n2[i]+1, n2[j]+1], q3mat[n3[i]+1, n3[j]+1]), k043)+Physics:-`*`(Physics:-`*`(q3mat[n2[i]+1, n2[j]+1], q4mat[n3[i]+1, n3[j]+1]), k034)+Physics:-`*`(Physics:-`*`(q2mat[n2[i]+1, n2[j]+1], q5mat[n3[i]+1, n3[j]+1]), k025)+Physics:-`*`(Physics:-`*`(q1mat[n2[i]+1, n2[j]+1], q6mat[n3[i]+1, n3[j]+1]), k016)), j = 1 .. nn), i = 1 .. nn)]); h_q23_m8 := Matrix(nn, nn, [seq(seq(Physics:-`*`(q0(n1[i], n1[j]), Physics:-`*`(Physics:-`*`(q7mat[n2[i]+1, n2[j]+1], q1mat[n3[i]+1, n3[j]+1]), k071)+Physics:-`*`(Physics:-`*`(q6mat[n2[i]+1, n2[j]+1], q2mat[n3[i]+1, n3[j]+1]), k062)+Physics:-`*`(Physics:-`*`(q5mat[n2[i]+1, n2[j]+1], q3mat[n3[i]+1, n3[j]+1]), k053)+Physics:-`*`(Physics:-`*`(q4mat[n2[i]+1, n2[j]+1], q4mat[n3[i]+1, n3[j]+1]), k044)+Physics:-`*`(Physics:-`*`(q3mat[n2[i]+1, n2[j]+1], q5mat[n3[i]+1, n3[j]+1]), k035)+Physics:-`*`(Physics:-`*`(q2mat[n2[i]+1, n2[j]+1], q6mat[n3[i]+1, n3[j]+1]), k026)+Physics:-`*`(Physics:-`*`(q1mat[n2[i]+1, n2[j]+1], q7mat[n3[i]+1, n3[j]+1]), k017)), j = 1 .. nn), i = 1 .. nn)]); h_q23 := `~`[Physics:-`*`](h_q23_m3+h_q23_m4+h_q23_m5+h_q23_m6+h_q23_m7+h_q23_m8, 219474.6); h_no := h_harm+h_q1+h_q2+h_q3; evals_no := linalg:-eigenvalues(h_no); v3no := evals_no[2]+Physics:-`*`(-1, evals_no[1]); v2no := evals_no[3]+Physics:-`*`(-1, evals_no[1]); v1no := evals_no[9]+Physics:-`*`(-1, evals_no[1]); uncoupled_freq := [v1no, v2no, v3no]; nq1 := kmax+1; nq2 := jmax+1; nq12 := Physics:-`*`(nq1, nq2); h_tot := h_harm+h_q1+h_q2+h_q3+h_q12+h_q13+h_q23; evals_pairmat := sort([linalg:-eigenvectors(h_tot)]); evals_pairlist := Vector(0); temp1 := 0; for j to 20 do temp1 := [seq(evals_pairmat[j, 3, 1][i], i = {1, 2, eval(nq1+1), eval(nq12+1)})]; seq(`if`(.50 <= abs(temp1[m]), Append(evals_pairlist, evals_pairmat[j, 1]), NULL), m = 1 .. nops(temp1)) end do; eval_size := ArrayTools:-Size(evals_pairlist, 1); eval_cont := [seq(evals_pairlist[i](2), i = 1 .. eval_size)]; if 4 < eval_size then ListTools:-Search(max(seq(eval_cont[i], i = 4 .. eval_size)), eval_cont); evals_pairlist := evals_pairlist[[1 .. 3, %]] end if; v3 := evals_pairlist[2](1)+Physics:-`*`(-1, evals_pairlist[1](1)); v2 := evals_pairlist[3](1)+Physics:-`*`(-1, evals_pairlist[1](1)); v1 := evals_pairlist[4](1)+Physics:-`*`(-1, evals_pairlist[1](1)); v3_t := evals_pairmat[2, 1]+Physics:-`*`(-1, evals_pairmat[1, 1]); v2_t := evals_pairmat[3, 1]+Physics:-`*`(-1, evals_pairmat[1, 1]); v1_t := evals_pairmat[9, 1]+Physics:-`*`(-1, evals_pairmat[1, 1]); coupled_freq := [v1, v2, v3]; coupled_freq_t := [v1_t, v2_t, v3_t]; return h_harm, h_q1, h_q2, h_q3, h_q12, h_q13, h_q23, uncoupled_freq, coupled_freq, coupled_freq_t, allqmat, row1, row2, row3, h_no, h_tot, evals_pairmat end proc

NULL

NULL

Download CoFr.mw

 

 

the attached, maybe(?)

  restart;
  interface(version);

`Standard Worksheet Interface, Maple 2019.0, Windows 7, March 9 2019 Build ID 1384062`

(1)

  with(LinearAlgebra):
  interface(rtablesize=20):
  with(plots):
  k:=2:M:=2:
  Tm:=(t,m)-> 2*t*Tm(t,m-1)-Tm(t,m-2);
  Tm(t,0):=1:
  Tm(t,1):=t:

  Tm2:=(tt,m)->subs(t=tt,Tm(t,m)):
  alpha:=(m)->piecewise(m=0,1/sqrt(Pi),sqrt(2)/sqrt(Pi));
  psii:=(n,m,x)->piecewise((n-1)/(2^(k-1)) <= x and x <= n/(2^(k-1)),
  alpha(m)*(2^(k/2))*Tm2(2^k*x-2*n+1,m), 0);
  psi:=(t)-> local i, j;
             Array([seq(seq(psii(i,j,t),j=0..M-1),i=1...2^(k-1))] ):
  for i from 1 to ((2^(k-1))*M) do
      r[i]:=evalf(psi((2*i-1)/((2^k)*M))):
  end do:
  m:=M*(2^(k-1)):
  xi:=(i,n)->((i+1)^(n+1)-2*i^(n+1)+(i-1)^(n+1));
  PB_n:= Matrix( m,
                 m,
                 (i, j)-> `if`( i=j,
                                1,
                                `if`( j>i,
                                      xi(j-1,n),
                                      0
                                    )
                               )
                );
  PB_n:=1/m^n/factorial(n+1)*PB_n;

Tm := proc (t, m) options operator, arrow; 2*t*Tm(t, m-1)-Tm(t, m-2) end proc

 

alpha := proc (m) options operator, arrow; piecewise(m = 0, 1/sqrt(Pi), sqrt(2)/sqrt(Pi)) end proc

 

psii := proc (n, m, x) options operator, arrow; piecewise((n-1)/2^(k-1) <= x and x <= n/2^(k-1), alpha(m)*2^((1/2)*k)*Tm2(2^k*x-2*n+1, m), 0) end proc

 

xi := proc (i, n) options operator, arrow; (i+1)^(n+1)-2*i^(n+1)+(i-1)^(n+1) end proc

 

Matrix(4, 4, {(1, 1) = 1, (1, 2) = 2^(n+1)-2, (1, 3) = 3^(n+1)-2*2^(n+1)+1, (1, 4) = 4^(n+1)-2*3^(n+1)+2^(n+1), (2, 1) = 0, (2, 2) = 1, (2, 3) = 3^(n+1)-2*2^(n+1)+1, (2, 4) = 4^(n+1)-2*3^(n+1)+2^(n+1), (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 4^(n+1)-2*3^(n+1)+2^(n+1), (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1})

 

Matrix(%id = 18446744074367853262)

(2)

 

Download getMat.mw

the following version is probably better

restart;
eq:=z=(r-4)/(r-2)/(r-6)+15:
sol:=solve(eq, r):
plot3d( sol[2], theta=0..2*Pi, z=0..40, coords=cylindrical, grid=[100,100], style=surface, axes=none);

 

 


 

Download SOR2.mw

upload a worksheet using the big green up-arrow in the Mapleprimes toolbar. No-one here feels like retyping your code: too time-consuming and error-prone.

I'd probably start with the command Student[Calculus1]:-SurfaceOfRevolution()

 

Solve() is an inert function - you probably(?) want solve() and maybe an allvalues() statement to evaluate the 'RootOfs', as in the attached


 

restart;
solve({(6*(-8*L*alpha^2*a[1]^3+2*L*a[1]^3-2*a[1]^3))*d-24*L*alpha^2*a[0]*a[1]^2+6*L*a[0]*a[1]^2-6*a[0]*a[1]^2+12*alpha^2*k^2*lambda*v*a[1]-12*alpha^2*k^2*lambda*a[1]-3*k^2*lambda*v*a[1]+3*k^2*lambda*a[1]+(3*(8*alpha^2*k^2*v^2*a[1]-16*alpha^2*k^2*v*a[1]+8*alpha^2*k^2*a[1]-2*k^2*v^2*a[1]+4*k^2*v*a[1]-2*k^2*a[1]))*d = 0, (-8*L*alpha^2*a[1]^3+2*L*a[1]^3-2*a[1]^3)*d^6+(-24*L*alpha^2*a[0]*a[1]^2+6*L*a[0]*a[1]^2-6*a[0]*a[1]^2)*d^5+(-24*L*alpha^2*a[0]^2*a[1]-24*L*alpha^2*a[1]^2*a[2]-4*alpha^4*a[1]+6*L*a[0]^2*a[1]+6*L*a[1]^2*a[2]-4*alpha^2*beta*a[1]+alpha^2*a[1]-6*a[0]^2*a[1]-6*a[1]^2*a[2]+beta*a[1])*d^4+(4*alpha^2*k^2*lambda*mu*a[1]-8*L*alpha^2*a[0]^3-48*L*alpha^2*a[0]*a[1]*a[2]-4*alpha^4*a[0]-k^2*lambda*mu*a[1]+2*L*a[0]^3+12*L*a[0]*a[1]*a[2]-4*alpha^2*beta*a[0]+alpha^2*a[0]-2*a[0]^3-12*a[0]*a[1]*a[2]+beta*a[0])*d^3+(-24*L*alpha^2*a[0]^2*a[2]-24*L*alpha^2*a[1]*a[2]^2-4*alpha^4*a[2]+6*L*a[0]^2*a[2]+6*L*a[1]*a[2]^2-4*alpha^2*beta*a[2]+alpha^2*a[2]-6*a[0]^2*a[2]-6*a[1]*a[2]^2+beta*a[2])*d^2+(-4*alpha^2*k^2*lambda*mu*a[2]-24*L*alpha^2*a[0]*a[2]^2+k^2*lambda*mu*a[2]+6*L*a[0]*a[2]^2-6*a[0]*a[2]^2)*d+8*alpha^2*k^2*mu^2*a[2]-2*a[2]^3-2*k^2*mu^2*a[2]-8*L*alpha^2*a[2]^3+2*L*a[2]^3 = 0, (6*(-8*L*alpha^2*a[1]^3+2*L*a[1]^3-2*a[1]^3))*d^5+(5*(-24*L*alpha^2*a[0]*a[1]^2+6*L*a[0]*a[1]^2-6*a[0]*a[1]^2))*d^4+(4*(-24*L*alpha^2*a[0]^2*a[1]-24*L*alpha^2*a[1]^2*a[2]-4*alpha^4*a[1]+6*L*a[0]^2*a[1]+6*L*a[1]^2*a[2]-4*alpha^2*beta*a[1]+alpha^2*a[1]-6*a[0]^2*a[1]-6*a[1]^2*a[2]+beta*a[1]))*d^3+(3*(4*alpha^2*k^2*lambda*mu*a[1]-8*L*alpha^2*a[0]^3-48*L*alpha^2*a[0]*a[1]*a[2]-4*alpha^4*a[0]-k^2*lambda*mu*a[1]+2*L*a[0]^3+12*L*a[0]*a[1]*a[2]-4*alpha^2*beta*a[0]+alpha^2*a[0]-2*a[0]^3-12*a[0]*a[1]*a[2]+beta*a[0]))*d^2+(4*alpha^2*k^2*lambda^2*a[1]+8*alpha^2*k^2*mu*v*a[1]-8*alpha^2*k^2*mu*a[1]-k^2*lambda^2*a[1]-2*k^2*mu*v*a[1]+2*k^2*mu*a[1])*d^3+(2*(-24*L*alpha^2*a[0]^2*a[2]-24*L*alpha^2*a[1]*a[2]^2-4*alpha^4*a[2]+6*L*a[0]^2*a[2]+6*L*a[1]*a[2]^2-4*alpha^2*beta*a[2]+alpha^2*a[2]-6*a[0]^2*a[2]-6*a[1]*a[2]^2+beta*a[2]))*d+12*alpha^2*k^2*lambda*mu*a[2]-24*L*alpha^2*a[0]*a[2]^2-3*k^2*lambda*mu*a[2]+6*L*a[0]*a[2]^2-6*a[0]*a[2]^2+(-4*alpha^2*k^2*lambda^2*a[2]-8*alpha^2*k^2*mu*v*a[2]+8*alpha^2*k^2*mu*a[2]+k^2*lambda^2*a[2]+2*k^2*mu*v*a[2]-2*k^2*mu*a[2])*d = 0, (15*(-8*L*alpha^2*a[1]^3+2*L*a[1]^3-2*a[1]^3))*d^2+(5*(-24*L*alpha^2*a[0]*a[1]^2+6*L*a[0]*a[1]^2-6*a[0]*a[1]^2))*d+(3*(12*alpha^2*k^2*lambda*v*a[1]-12*alpha^2*k^2*lambda*a[1]-3*k^2*lambda*v*a[1]+3*k^2*lambda*a[1]))*d+(3*(8*alpha^2*k^2*v^2*a[1]-16*alpha^2*k^2*v*a[1]+8*alpha^2*k^2*a[1]-2*k^2*v^2*a[1]+4*k^2*v*a[1]-2*k^2*a[1]))*d^2+8*alpha^2*k^2*mu*v*a[1]-6*a[0]^2*a[1]-6*a[1]^2*a[2]-4*alpha^4*a[1]+4*alpha^2*k^2*lambda^2*a[1]-8*alpha^2*k^2*mu*a[1]-2*k^2*mu*v*a[1]-24*L*alpha^2*a[0]^2*a[1]-24*L*alpha^2*a[1]^2*a[2]-k^2*lambda^2*a[1]+2*k^2*mu*a[1]+6*L*a[0]^2*a[1]+6*L*a[1]^2*a[2]-4*alpha^2*beta*a[1]+alpha^2*a[1]+beta*a[1] = 0, (15*(-8*L*alpha^2*a[1]^3+2*L*a[1]^3-2*a[1]^3))*d^4+(10*(-24*L*alpha^2*a[0]*a[1]^2+6*L*a[0]*a[1]^2-6*a[0]*a[1]^2))*d^3+(6*(-24*L*alpha^2*a[0]^2*a[1]-24*L*alpha^2*a[1]^2*a[2]-4*alpha^4*a[1]+6*L*a[0]^2*a[1]+6*L*a[1]^2*a[2]-4*alpha^2*beta*a[1]+alpha^2*a[1]-6*a[0]^2*a[1]-6*a[1]^2*a[2]+beta*a[1]))*d^2+(12*alpha^2*k^2*lambda*v*a[1]-12*alpha^2*k^2*lambda*a[1]-3*k^2*lambda*v*a[1]+3*k^2*lambda*a[1])*d^3+(3*(4*alpha^2*k^2*lambda*mu*a[1]-8*L*alpha^2*a[0]^3-48*L*alpha^2*a[0]*a[1]*a[2]-4*alpha^4*a[0]-k^2*lambda*mu*a[1]+2*L*a[0]^3+12*L*a[0]*a[1]*a[2]-4*alpha^2*beta*a[0]+alpha^2*a[0]-2*a[0]^3-12*a[0]*a[1]*a[2]+beta*a[0]))*d+(3*(4*alpha^2*k^2*lambda^2*a[1]+8*alpha^2*k^2*mu*v*a[1]-8*alpha^2*k^2*mu*a[1]-k^2*lambda^2*a[1]-2*k^2*mu*v*a[1]+2*k^2*mu*a[1]))*d^2+(-12*alpha^2*k^2*lambda*v*a[2]+12*alpha^2*k^2*lambda*a[2]+3*k^2*lambda*v*a[2]-3*k^2*lambda*a[2])*d+8*alpha^2*k^2*mu*v*a[2]-6*a[0]^2*a[2]-6*a[1]*a[2]^2-4*alpha^4*a[2]+4*alpha^2*k^2*lambda^2*a[2]-8*alpha^2*k^2*mu*a[2]-2*k^2*mu*v*a[2]-24*L*alpha^2*a[0]^2*a[2]-24*L*alpha^2*a[1]*a[2]^2-k^2*lambda^2*a[2]+2*k^2*mu*a[2]+6*L*a[0]^2*a[2]+6*L*a[1]*a[2]^2-4*alpha^2*beta*a[2]+alpha^2*a[2]+beta*a[2] = 0, (20*(-8*L*alpha^2*a[1]^3+2*L*a[1]^3-2*a[1]^3))*d^3+(10*(-24*L*alpha^2*a[0]*a[1]^2+6*L*a[0]*a[1]^2-6*a[0]*a[1]^2))*d^2+(4*(-24*L*alpha^2*a[0]^2*a[1]-24*L*alpha^2*a[1]^2*a[2]-4*alpha^4*a[1]+6*L*a[0]^2*a[1]+6*L*a[1]^2*a[2]-4*alpha^2*beta*a[1]+alpha^2*a[1]-6*a[0]^2*a[1]-6*a[1]^2*a[2]+beta*a[1]))*d+(8*alpha^2*k^2*v^2*a[1]-16*alpha^2*k^2*v*a[1]+8*alpha^2*k^2*a[1]-2*k^2*v^2*a[1]+4*k^2*v*a[1]-2*k^2*a[1])*d^3+(3*(4*alpha^2*k^2*lambda^2*a[1]+8*alpha^2*k^2*mu*v*a[1]-8*alpha^2*k^2*mu*a[1]-k^2*lambda^2*a[1]-2*k^2*mu*v*a[1]+2*k^2*mu*a[1]))*d+(3*(12*alpha^2*k^2*lambda*v*a[1]-12*alpha^2*k^2*lambda*a[1]-3*k^2*lambda*v*a[1]+3*k^2*lambda*a[1]))*d^2+(-8*alpha^2*k^2*v^2*a[2]+16*alpha^2*k^2*v*a[2]-8*alpha^2*k^2*a[2]+2*k^2*v^2*a[2]-4*k^2*v*a[2]+2*k^2*a[2])*d+4*alpha^2*k^2*lambda*v*a[2]+2*L*a[0]^3-4*alpha^4*a[0]-48*L*alpha^2*a[0]*a[1]*a[2]+4*alpha^2*k^2*lambda*mu*a[1]-4*alpha^2*k^2*lambda*a[2]-k^2*lambda*v*a[2]-k^2*lambda*mu*a[1]+12*L*a[0]*a[1]*a[2]+k^2*lambda*a[2]-8*L*alpha^2*a[0]^3-12*a[0]*a[1]*a[2]-4*alpha^2*beta*a[0]+alpha^2*a[0]+beta*a[0]-2*a[0]^3 = 0, 8*alpha^2*k^2*v^2*a[1]-8*L*alpha^2*a[1]^3-16*alpha^2*k^2*v*a[1]+8*alpha^2*k^2*a[1]-2*k^2*v^2*a[1]+2*L*a[1]^3+4*k^2*v*a[1]-2*k^2*a[1]-2*a[1]^3 = 0}, {alpha, a[0], a[1], a[2]});

{alpha = alpha, a[0] = 0, a[1] = 0, a[2] = 0}, {alpha = alpha, a[0] = RootOf((8*L*alpha^2-2*L+2)*_Z^2+4*alpha^4+4*alpha^2*beta-alpha^2-beta), a[1] = 0, a[2] = 0}, {alpha = RootOf(k^2*lambda^2-4*k^2*mu*v+4*k^2*mu+2*_Z^2+2*beta), a[0] = RootOf((8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4)*_Z^2-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)*(2*d*v-2*d-lambda)*k, a[1] = 0, a[2] = -2*(d^2*v-d^2-d*lambda+mu)*RootOf((8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4)*_Z^2-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)*k}, {alpha = RootOf(k^2*lambda^2-4*k^2*mu*v+4*k^2*mu+2*_Z^2+2*beta), a[0] = -(1/2)*k*(4*d*k^2*lambda^2*v-16*d*k^2*mu*v^2-4*d*k^2*lambda^2+32*d*k^2*mu*v-2*k^2*lambda^3+8*k^2*lambda*mu*v-16*d*k^2*mu-8*k^2*lambda*mu+8*beta*d*v-8*beta*d-4*beta*lambda+2*d*v-2*d-lambda)/(RootOf(_Z^2*(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1)-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)*(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1)), a[1] = RootOf(_Z^2*(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1)-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)*(v-1)*k, a[2] = 0}

(1)

allvalues~([%]);

[{alpha = alpha, a[0] = 0, a[1] = 0, a[2] = 0}, {alpha = alpha, a[0] = (-(4*alpha^4+4*alpha^2*beta-alpha^2-beta)/(8*L*alpha^2-2*L+2))^(1/2), a[1] = 0, a[2] = 0}, {alpha = alpha, a[0] = -(-(4*alpha^4+4*alpha^2*beta-alpha^2-beta)/(8*L*alpha^2-2*L+2))^(1/2), a[1] = 0, a[2] = 0}, {alpha = (-(1/2)*k^2*lambda^2+2*k^2*mu*v-2*k^2*mu-beta)^(1/2), a[0] = (-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4))^(1/2)*(2*d*v-2*d-lambda)*k, a[1] = 0, a[2] = -2*(d^2*v-d^2-d*lambda+mu)*(-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4))^(1/2)*k}, {alpha = (-(1/2)*k^2*lambda^2+2*k^2*mu*v-2*k^2*mu-beta)^(1/2), a[0] = -(-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4))^(1/2)*(2*d*v-2*d-lambda)*k, a[1] = 0, a[2] = 2*(d^2*v-d^2-d*lambda+mu)*(-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4))^(1/2)*k}, {alpha = -(-(1/2)*k^2*lambda^2+2*k^2*mu*v-2*k^2*mu-beta)^(1/2), a[0] = (-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4))^(1/2)*(2*d*v-2*d-lambda)*k, a[1] = 0, a[2] = -2*(d^2*v-d^2-d*lambda+mu)*(-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4))^(1/2)*k}, {alpha = -(-(1/2)*k^2*lambda^2+2*k^2*mu*v-2*k^2*mu-beta)^(1/2), a[0] = -(-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4))^(1/2)*(2*d*v-2*d-lambda)*k, a[1] = 0, a[2] = 2*(d^2*v-d^2-d*lambda+mu)*(-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(8*L*k^2*lambda^2-32*L*k^2*mu*v+32*L*k^2*mu+16*L*beta+4*L-4))^(1/2)*k}, {alpha = (-(1/2)*k^2*lambda^2+2*k^2*mu*v-2*k^2*mu-beta)^(1/2), a[0] = -(1/2)*k*(4*d*k^2*lambda^2*v-16*d*k^2*mu*v^2-4*d*k^2*lambda^2+32*d*k^2*mu*v-2*k^2*lambda^3+8*k^2*lambda*mu*v-16*d*k^2*mu-8*k^2*lambda*mu+8*beta*d*v-8*beta*d-4*beta*lambda+2*d*v-2*d-lambda)/((-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1))^(1/2)*(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1)), a[1] = (-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1))^(1/2)*(v-1)*k, a[2] = 0}, {alpha = (-(1/2)*k^2*lambda^2+2*k^2*mu*v-2*k^2*mu-beta)^(1/2), a[0] = (1/2)*k*(4*d*k^2*lambda^2*v-16*d*k^2*mu*v^2-4*d*k^2*lambda^2+32*d*k^2*mu*v-2*k^2*lambda^3+8*k^2*lambda*mu*v-16*d*k^2*mu-8*k^2*lambda*mu+8*beta*d*v-8*beta*d-4*beta*lambda+2*d*v-2*d-lambda)/((-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1))^(1/2)*(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1)), a[1] = -(-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1))^(1/2)*(v-1)*k, a[2] = 0}, {alpha = -(-(1/2)*k^2*lambda^2+2*k^2*mu*v-2*k^2*mu-beta)^(1/2), a[0] = -(1/2)*k*(4*d*k^2*lambda^2*v-16*d*k^2*mu*v^2-4*d*k^2*lambda^2+32*d*k^2*mu*v-2*k^2*lambda^3+8*k^2*lambda*mu*v-16*d*k^2*mu-8*k^2*lambda*mu+8*beta*d*v-8*beta*d-4*beta*lambda+2*d*v-2*d-lambda)/((-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1))^(1/2)*(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1)), a[1] = (-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1))^(1/2)*(v-1)*k, a[2] = 0}, {alpha = -(-(1/2)*k^2*lambda^2+2*k^2*mu*v-2*k^2*mu-beta)^(1/2), a[0] = (1/2)*k*(4*d*k^2*lambda^2*v-16*d*k^2*mu*v^2-4*d*k^2*lambda^2+32*d*k^2*mu*v-2*k^2*lambda^3+8*k^2*lambda*mu*v-16*d*k^2*mu-8*k^2*lambda*mu+8*beta*d*v-8*beta*d-4*beta*lambda+2*d*v-2*d-lambda)/((-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1))^(1/2)*(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1)), a[1] = -(-(-2*k^2*lambda^2+8*k^2*mu*v-8*k^2*mu-4*beta-1)/(2*L*k^2*lambda^2-8*L*k^2*mu*v+8*L*k^2*mu+4*L*beta+L-1))^(1/2)*(v-1)*k, a[2] = 0}]

(2)

 


 

Download doSolve.mw

 

Maple's plotting routines are adaptive.

There is no guarantee that the number of points used to generate different curves will be the same

Even if number of points is the same, there is no guarantee that the same x-values have been used to produce all the plots. So you cannot take the x-values from a plot of (x1, y1) and uses thes to generate a matrix with columns x1, y1, y2, y3, y4, y5, y6, y7, where y2, y3, y4, y5, y6, y7 arise from different plots

You have two choices

  1. If you really want to extract data from several Maple plots and export them, you will need to do this into several different data files, each of which contains its own (unique) x-y data. Using this method you will get seven output files, each of which contains (about) 200 rows and 2 columns
  2. On the other hand, you could generate an output matrix directly (ie without extracting data from plots) simply by calculating the values for each function at the same set of x-values. In this way one can get a single matrix, with x-values in the first column and function values in the subsequent seven columns

The attached implements both of the above approaches  in the final couple of execution groups.

You will have to change the filepath used in the attached to something appropriate for your machine


 

restart

u := -.628767281517651653319062551564*t-0.668228578765714691262317435292e-2*t^2+.214453812334969261759874843315*t^3-0.754482317038267703413764826760e-1*t^4+.101654092326534854140522092073*t^5-.230150346728028966545611699989*t^6+.280322383361924053601807027731*t^7-.247903345761556356977695694547*t^8+.178410936157224617321140075295*t^9-.103903421341449622624578629693*t^10+.590877216504511906626000249484+0.475644216784992732878138002113e-1*t^11-0.168193324621169957791081606619e-1*t^12+0.455321539911051144632769231901e-2*t^13-0.934897021209409120472374464828e-3*t^14+0.143227945580893189263690052053e-3*t^15-0.158869689089997293849868657822e-4*t^16+0.120728607755254184751628676867e-5*t^17-5.63137388794357080952813577816*10^(-8)*t^18+1.21736603352619990127647551648*10^(-9)*t^19; U := eval(u, t = t^alpha/GAMMA(alpha+1))

-.628767281517651653319062551564*t-0.668228578765714691262317435292e-2*t^2+.214453812334969261759874843315*t^3-0.754482317038267703413764826760e-1*t^4+.101654092326534854140522092073*t^5-.230150346728028966545611699989*t^6+.280322383361924053601807027731*t^7-.247903345761556356977695694547*t^8+.178410936157224617321140075295*t^9-.103903421341449622624578629693*t^10+.590877216504511906626000249484+0.475644216784992732878138002113e-1*t^11-0.168193324621169957791081606619e-1*t^12+0.455321539911051144632769231901e-2*t^13-0.934897021209409120472374464828e-3*t^14+0.143227945580893189263690052053e-3*t^15-0.158869689089997293849868657822e-4*t^16+0.120728607755254184751628676867e-5*t^17-0.5631373888e-7*t^18+0.1217366034e-8*t^19

 

-.628767281517651653319062551564*t^alpha/GAMMA(alpha+1)-0.668228578765714691262317435292e-2*(t^alpha)^2/GAMMA(alpha+1)^2+.214453812334969261759874843315*(t^alpha)^3/GAMMA(alpha+1)^3-0.754482317038267703413764826760e-1*(t^alpha)^4/GAMMA(alpha+1)^4+.101654092326534854140522092073*(t^alpha)^5/GAMMA(alpha+1)^5-.230150346728028966545611699989*(t^alpha)^6/GAMMA(alpha+1)^6+.280322383361924053601807027731*(t^alpha)^7/GAMMA(alpha+1)^7-.247903345761556356977695694547*(t^alpha)^8/GAMMA(alpha+1)^8+.178410936157224617321140075295*(t^alpha)^9/GAMMA(alpha+1)^9-.103903421341449622624578629693*(t^alpha)^10/GAMMA(alpha+1)^10+.590877216504511906626000249484+0.475644216784992732878138002113e-1*(t^alpha)^11/GAMMA(alpha+1)^11-0.168193324621169957791081606619e-1*(t^alpha)^12/GAMMA(alpha+1)^12+0.455321539911051144632769231901e-2*(t^alpha)^13/GAMMA(alpha+1)^13-0.934897021209409120472374464828e-3*(t^alpha)^14/GAMMA(alpha+1)^14+0.143227945580893189263690052053e-3*(t^alpha)^15/GAMMA(alpha+1)^15-0.158869689089997293849868657822e-4*(t^alpha)^16/GAMMA(alpha+1)^16+0.120728607755254184751628676867e-5*(t^alpha)^17/GAMMA(alpha+1)^17-0.5631373888e-7*(t^alpha)^18/GAMMA(alpha+1)^18+0.1217366034e-8*(t^alpha)^19/GAMMA(alpha+1)^19

(1)

T := -1.59721123452136569934800873143*t+.140914158578234612583524954189*t^2+.184326725551320973168311890071*t^3+3.27679250886063610552711558733*t^4-9.19121403283618469299869253913*t^5+16.2210977323449252708580078591*t^6-22.6274998117906738354944873304*t^7+24.6405998730600422897902054806*t^8-20.4642695873215256847411520249*t^9+12.9399237576588503716958279801*t^10-6.27095040567535707739636899345*t^11+2.34164974920843825133818318661*t^12-.673785216636421465182335202225*t^13+.148302898808793236284492098988*t^14-0.245445239997071611634751192522e-1*t^15+0.295919502140931594969136594280e-2*t^16-0.245522039023190415694499765818e-3*t^17+0.125414172181754774790292035344e-4*t^18-2.97429482868330719071713325159*10^(-7)*t^19+1.00000000000000000000000000001; TT := eval(T, t = t^alpha/GAMMA(alpha+1))

-1.59721123452136569934800873143*t+.140914158578234612583524954189*t^2+.184326725551320973168311890071*t^3+3.27679250886063610552711558733*t^4-9.19121403283618469299869253913*t^5+16.2210977323449252708580078591*t^6-22.6274998117906738354944873304*t^7+24.6405998730600422897902054806*t^8-20.4642695873215256847411520249*t^9+12.9399237576588503716958279801*t^10-6.27095040567535707739636899345*t^11+2.34164974920843825133818318661*t^12-.673785216636421465182335202225*t^13+.148302898808793236284492098988*t^14-0.245445239997071611634751192522e-1*t^15+0.295919502140931594969136594280e-2*t^16-0.245522039023190415694499765818e-3*t^17+0.125414172181754774790292035344e-4*t^18-0.2974294829e-6*t^19+1.00000000000000000000000000001

 

-1.59721123452136569934800873143*t^alpha/GAMMA(alpha+1)+.140914158578234612583524954189*(t^alpha)^2/GAMMA(alpha+1)^2+.184326725551320973168311890071*(t^alpha)^3/GAMMA(alpha+1)^3+3.27679250886063610552711558733*(t^alpha)^4/GAMMA(alpha+1)^4-9.19121403283618469299869253913*(t^alpha)^5/GAMMA(alpha+1)^5+16.2210977323449252708580078591*(t^alpha)^6/GAMMA(alpha+1)^6-22.6274998117906738354944873304*(t^alpha)^7/GAMMA(alpha+1)^7+24.6405998730600422897902054806*(t^alpha)^8/GAMMA(alpha+1)^8-20.4642695873215256847411520249*(t^alpha)^9/GAMMA(alpha+1)^9+12.9399237576588503716958279801*(t^alpha)^10/GAMMA(alpha+1)^10-6.27095040567535707739636899345*(t^alpha)^11/GAMMA(alpha+1)^11+2.34164974920843825133818318661*(t^alpha)^12/GAMMA(alpha+1)^12-.673785216636421465182335202225*(t^alpha)^13/GAMMA(alpha+1)^13+.148302898808793236284492098988*(t^alpha)^14/GAMMA(alpha+1)^14-0.245445239997071611634751192522e-1*(t^alpha)^15/GAMMA(alpha+1)^15+0.295919502140931594969136594280e-2*(t^alpha)^16/GAMMA(alpha+1)^16-0.245522039023190415694499765818e-3*(t^alpha)^17/GAMMA(alpha+1)^17+0.125414172181754774790292035344e-4*(t^alpha)^18/GAMMA(alpha+1)^18-0.2974294829e-6*(t^alpha)^19/GAMMA(alpha+1)^19+1.00000000000000000000000000001

(2)

U1 := eval(U, alpha = .5); U2 := eval(U, alpha = .6); U3 := eval(U, alpha = .7); U4 := eval(U, alpha = .8); U5 := eval(U, alpha = .9); U6 := eval(U, alpha = 1); U77 := 100000000000000000000000000000000000000000000*erf((1/250000000000000)*sqrt(891963910829259026055675914205))*sqrt(Pi)/(100000000000000000000000000000000000000000000*erf((1/250000000000000)*sqrt(891963910829259026055675914205))*sqrt(Pi)+130134883134501207078143767261*sqrt(891963910829259026055675914205))-100000000000000000000000000000000000000000000*erf((1/1000000000000000)*sqrt(891963910829259026055675914205)*eta)*sqrt(Pi)/(100000000000000000000000000000000000000000000*erf((1/250000000000000)*sqrt(891963910829259026055675914205))*sqrt(Pi)+130134883134501207078143767261*sqrt(891963910829259026055675914205)); U7 := eval(U77, eta = t); T1 := eval(TT, alpha = .5); T2 := eval(TT, alpha = .6); T3 := eval(TT, alpha = .7); T4 := eval(TT, alpha = .8); T5 := eval(TT, alpha = .9); T6 := eval(TT, alpha = 1); T77 := 1-erf((1/8555318768232570228568063446)*sqrt(144545018049207946484247493468681726749763693119033260670)*eta)/erf((2/4277659384116285114284031723)*sqrt(144545018049207946484247493468681726749763693119033260670)); T7 := eval(T77, eta = t)

100000000000000000000000000000000000000000000*erf((1/250000000000000)*891963910829259026055675914205^(1/2))*Pi^(1/2)/(100000000000000000000000000000000000000000000*erf((1/250000000000000)*891963910829259026055675914205^(1/2))*Pi^(1/2)+130134883134501207078143767261*891963910829259026055675914205^(1/2))-100000000000000000000000000000000000000000000*erf((1/1000000000000000)*891963910829259026055675914205^(1/2)*t)*Pi^(1/2)/(100000000000000000000000000000000000000000000*erf((1/250000000000000)*891963910829259026055675914205^(1/2))*Pi^(1/2)+130134883134501207078143767261*891963910829259026055675914205^(1/2))

 

1-erf((1/8555318768232570228568063446)*144545018049207946484247493468681726749763693119033260670^(1/2)*t)/erf((2/4277659384116285114284031723)*144545018049207946484247493468681726749763693119033260670^(1/2))

(3)

plot([U1, U2, U3, U4, U5, U6, U7], t = 0 .. 4, color = [red, blue, green, cyan, black, purple, "DeepPink"]); S1 := plot([U1, U2, U3, U4, U5, U6, U7], t = 0 .. 4, color = [red, blue, green, cyan, black, purple, "DeepPink"]); plot([T1, T2, T3, T4, T5, T6, T7], t = 0 .. 4, color = [red, blue, green, cyan, black, purple, "DeepPink"]); S2 := plot([T1, T2, T3, T4, T5, T6, T7], t = 0 .. 4, color = [red, blue, T1green, cyan, black, purple, "DeepPink"])

 

 

with(plots):

 

p11 := display({S2}, axes = boxed, thickness = 2); dat1 := `~`[plottools:-getdata]([p11]); for i while i <= 7 do A[i] := dat1[i, 3]; Y[i] := A[i][1 .. -1, 2] end do; X := A[1][1 .. -1, 1]; MM1 := `<|>`(X, `$`(Y[j2], j2 = 1 .. 7))

 

[["curve", [0. .. 4., -3.13342742419564502*10^(-7) .. 1.], _rtable[18446744074580050518]], ["curve", [0. .. 4., -0.655056087701666456e-4 .. 1.], _rtable[18446744074580050638]], ["curve", [0. .. 4., 0.884660748723398171e-5 .. 1.], _rtable[18446744074580050758]], ["curve", [0. .. 4., 0.424814094387437891e-4 .. 1.], _rtable[18446744074580050878]], ["curve", [0. .. 4., -0.436799487579264678e-2 .. 1.], _rtable[18446744074580050998]], ["curve", [0. .. 4., -0.130796688608825207e-4 .. 1.], _rtable[18446744074580051118]], ["curve", [0. .. 4., 0. .. 1.], _rtable[18446744074580051238]]]

 

_rtable[18446744074580050518]

 

_rtable[18446744074580022926]

 

_rtable[18446744074580050638]

 

_rtable[18446744074580019910]

 

_rtable[18446744074580050758]

 

_rtable[18446744074580016894]

 

_rtable[18446744074580050878]

 

_rtable[18446744074580005702]

 

_rtable[18446744074580050998]

 

_rtable[18446744074580002686]

 

_rtable[18446744074580051118]

 

_rtable[18446744074579999670]

 

_rtable[18446744074580051238]

 

_rtable[18446744074424885246]

 

Error, (in Matrix) this entry is too tall or too short: Vector(204, {(1) = 1.0, (2) = .9236715124867816, (3) = .8846539212704589, (4) = .853296723180952, (5) = .8261372152465285, (6) = .7850982395717462, (7) = .7493208608333364, (8) = .6809566063078758, (9) = .6232461324642893, (10) = .5733608398336795, (11) = .5322087236589063, (12) = .4937540687102193, (13) = .457693161090617, (14) = .4249445998228353, (15) = .3941405636373234, (16) = .3691419589560857, (17) = .3431280736389174, (18) = .3190444867517285, (19) = .2975764708459453, (20) = .2794286709540693, (21) = .25937298220622207, (22) = .24357288561927237, (23) = .22645 ... 799450080739e-4, (204) = -0.6550560877016665e-4}, datatype = float[8])

 

ExportMatrix("C:/Users/DELL/Dropbox///SP alpha varies u.dat", MM1)

Error, invalid input: ExportMatrix expects its 2nd argument, M, to be of type {Matrix, list(Matrix)}, but received MM1

 

 

_rtable[18446744074424763446]

 

Matrix(%id = 18446744074424758990)

(4)

#
# Return the dimensions of each curve matrix, just to
# check.
#
# Note that these are different!
#
# That's adaptive plotting routines for you!
#
  seq( [op(1, plottools:-getdata(S2)[j][3])], j=1..7);
#
# So the following will output 7 dataFiles, each
# of which contains x-y data
#
  seq( ExportMatrix
       ( cat( "C:/Users/TomLeslie/Desktop/datFile", j, ".dat"),
         plottools:-getdata(S2)[j][3]
       ),
       j=1..7
     );

[206, 2], [204, 2], [202, 2], [200, 2], [200, 2], [200, 2], [200, 2]

 

9889, 9843, 9696, 9600, 9688, 9688, 9600

(5)

#
# If the objective is to produce a 200 X 8 matrix
# with columns x, T1(x)......T7(x), where x-values
# in column1 are equally spaced, and range from
# 0..1, then this can be achieved directly with the
# following.
#
  T[1] := eval(TT, alpha = 0.5):
  T[2] := eval(TT, alpha = 0.6):
  T[3] := eval(TT, alpha = 0.7):
  T[4] := eval(TT, alpha = 0.8):
  T[5] := eval(TT, alpha = 0.9):
  T[6] := eval(TT, alpha = 1):
  T77 := 1 - erf(1/8555318768232570228568063446*sqrt(144545018049207946484247493468681726749763693119033260670)*eta)/erf(2/4277659384116285114284031723*sqrt(144545018049207946484247493468681726749763693119033260670)):
  T[7] := eval(T77, eta = t):
  MM:= Matrix
       ( npts,
         8,
         (i, j)->`if`( j=1,
                       evalf((i-1)/(npts-1)),
                       evalf( eval(T[j-1], t=(i-1)/(npts-1)))
                     )
      );
#
# Matrix can then be output directly
#
  ExportMatrix( "C:/Users/TomLeslie/Desktop/AllDataFile.dat", MM);
#
# If necessary one can produce plots directly from this matrix
#
  cs := [red, blue, green, cyan, black, purple, "DeepPink"]:
  display( [seq( pointplot( MM[..,1], MM[..,j+1], color=cs[j]), j=1..7)] );

_rtable[18446744074424777054]

 

20048

 

 

 

 

 

 


 

Download saveMatrix.mw

what you mean by a "complete trigonometric circle".

However the attached replicates the image you have posted. This may (or may not!) be useful.

BTW - this plot renders somewhat better in a Maple worksheet, than it does on this site

  restart;
  with(plots):
  doSQ:= proc( r, ang, col )
               uses plottools:
               local i, j,
                     x1:=-r*cos(ang),
                     y1:=-r*sin(ang),
                     P:= Array( 1..4,
                                [ [  x1,  y1],
                                  [  x1, -y1],
                                  [ -x1,  y1],
                                  [ -x1, -y1]
                                ]
                              ):
               return seq
                      ( seq
                        ( line
                          ( P[i],
                            P[j],
                            color=col
                          ),
                          j=i+1..4
                        ),
                        i=1..4
                      );
        end proc:
  display( [ polarplot(1),
             doSQ(1, Pi/6, red),
             doSQ(1, Pi/4, green),
             doSQ(1, Pi/3, blue)
           ],
           size=[800,800],
           scaling=constrained
         );

 

 

NULL


 

Download trigCir.mw

the attached.

Note that (as Carl has pointed out), the firts two operations yield the same result

  EQ := V^3-R*T*V^2/P-(B^2+P*B*R*T*sqrt(T)/(P-A))*V-A*B/(P*sqrt(T)) = 0;
  diff(EQ,T,V);
  diff(EQ,V,T);
  diff(EQ,V,V);

V^3-R*T*V^2/P-(B^2+P*B*R*T^(3/2)/(P-A))*V-A*B/(P*T^(1/2)) = 0

 

-2*R*V/P-(3/2)*P*B*R*T^(1/2)/(P-A) = 0

 

-2*R*V/P-(3/2)*P*B*R*T^(1/2)/(P-A) = 0

 

6*V-2*R*T/P = 0

(1)

 

 


 

Download doDiffs.mw

Command has been "repositioned".

It is now under the "Evaluate" tab

No idea why.

You have defined T[1], T[2] etc as functions. If you wish to evaluate these for some argument, then you need to use T[1](argumentValue). For example

T[1](3);
T[2](7);
or even
T[1](z)

Please stop posting "pictures" of code. Post editable/executable code by uploading your worksheet using the big green up-arrow in the Mapleprimes toolbar

There are several syntax errors, which essentially prevent your worksheet from executing. I have fixed these in the attached.

However (I am nearly certain that) there are further logical errors, which means that this code will not produce what you wish. You need to consider

  1. Do you expect the local variable 'n' in the procedure and the variable 'n' used at the top-level to be the same entity - because they are not! If you wish to evaluate the variable sigma_xx at the value n=loop_index, then the simplest way is probably to change the loop index to 'j' and write eval( sigma_xx, n=j). I have made this change in the second execution group of the attached.
  2. You have two for loops - both of these assign something to Fx, Should one of these be Fy? I have made this change in the second execution group of the attached
  3. Each time either of these for loops is executed, it will overwrite the value obtained on the previous iteration. I have adjusted the code so that  each for loop generates a sequence of values for Fx and Fy
  4. the return statement needs to be modified to return both sequences.

With these changes (some/all of which may not be what you actually want, becuase I'm guessing!), the code executes, and produces values for some of the integrals, although some are 'undefined'. I don't propose to work out why the integration sometimes return undefined values until/unless you confirm the validity of the changes I have made

restart;
F0:=proc(sigma__xx,N)
local  x,y,Fx,Fy: #removed superfluous comma
assume (w,real,w>0):assume (h,real,h>0):

for n from 0 by 2 to N do Fx:=integrate(sigma__xx*cos(w*Zeta*h),Zeta=0..infinity):
end do;
for n from 1 by 2 to N do Fx:=integrate(sigma__xx*sin(w*Zeta*h),Zeta=0..infinity):
end do;
return [Fx]:
end proc;

sigma__xx := -(sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2+sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n-2*sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h)))*Zeta^2*h^2/(Zeta^2*h^2+h^2)^2+(sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2-3*sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n); #added closing parenthesis (may be in wrong place!)

F0(sigma__xx,4); # 'for' loops in proc need numeric final value


 

Warning, `n` is implicitly declared local to procedure `F0`

 

proc (sigma__xx, N) local x, y, Fx, Fy, n; assume(w, real, 0 < w); assume(h, real, 0 < h); for n from 0 by 2 to N do Fx := integrate(sigma__xx*cos(w*Zeta*h), Zeta = 0 .. infinity) end do; for n by 2 to N do Fx := integrate(sigma__xx*sin(w*Zeta*h), Zeta = 0 .. infinity) end do; return [Fx] end proc

 

-(((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2+((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h))*n-2*((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h)))*Zeta^2*h^2/(Zeta^2*h^2+h^2)^2+((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2-3*((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(h, Zeta*h))*n

 

[int((-(((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta))*n^2+((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta))*n-2*((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta)))*Zeta^2*h^2/(Zeta^2*h^2+h^2)^2+((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta))*n^2-3*((Zeta^2*h^2+h^2)^(1/2))^(-n+2)*cos(n*arctan(1, Zeta))*n)*sin(w*Zeta*h), Zeta = 0 .. infinity)]

(1)

restart;
F0:= proc(sigma__xx,N)
          local  j, x, y, Fx:=NULL, Fy:=NULL:
          assume (w,real,w>0):
          assume (h,real,h>0):
          for j from 0 by 2 to N do
              Fx:= Fx,
                   integrate
                   ( eval
                     ( sigma__xx,
                       n=j
                     )*cos(w*Zeta*h),
                     Zeta = 0..infinity
                   ):
          end do;
          for j from 1 by 2 to N do
              Fy:= Fy,
                   integrate
                   ( eval
                     ( sigma__xx,
                       n=j
                     )*sin(w*Zeta*h),
                     Zeta = 0..infinity
                   ):
          end do;
          return [Fx], [Fy]:
     end proc:

sigma__xx := -(sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2+sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n-2*sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h)))*Zeta^2*h^2/(Zeta^2*h^2+h^2)^2+(sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n^2-3*sqrt(Zeta^2*h^2+h^2)^(-n+2)*cos(n*arctan(h, Zeta*h))*n):

F0(sigma__xx,5);

[undefined, undefined, (1/16)*exp(-w*h)*Pi*(3*h^4*w^4+32*h^4*w^2-18*h^3*w^3-32*h^3*w+9*h^2*w^2+9*h*w+9)/h^4], [undefined, -(5/24)*w*exp(-w*h)*Pi*(2*h^2*w^2-9*h*w+3)/h^2, (1/240)*w*exp(-w*h)*Pi*(14*h^4*w^4+400*h^4*w^2-105*h^3*w^3-600*h^3*w+70*h^2*w^2+105*h*w+105)/h^4]

(2)

 

 


 

Download doInt.mw

is that you want to compare the results of solving the ODE using Maple's numeric solver with those obtained from your numeric process. And compare both of these with Maple's exact solution.

I have modified your code to do this in the attached. The Maple numeric solver (on default accuracy settings) is slightly more accurate than your numeric process - but not by much!

restart;
#Differential Equation: 
#y''(t)+a(t)y'+b(t)y=f(t)
#IC's:
#y(0)=alpha,y'(0)=beta

k:=2:  
M:=4:  
N:=2^(k-1)*M:
#  Let's define Legendre Polynomials
P[0]:= t -> 1:
P[1]:= t -> t:
for m from 2 to M-1 do
   P[m]:= unapply(expand((2*m-1)*t*P[m-1](t)/m - (m-1)*P[m-2](t)/m), t)
end do:

 
L:=proc(n,m) local a,b;
  a := (2*n-2)/2^k;
  b := 2*n/2^k;
  return piecewise(a <= t and t <= b, P[m](2^k*t-2*n+1)*sqrt(((m+1)/2)*2^k))
end proc:

Ld := Vector(N):
LL := Matrix(N, N):
T := Vector(N):
for j from 1 to N do T[j] := (j-1/2)/N end do:
for n from 1 to 2^(k-1) do
  for m from 0 to M-1 do
i := n+2^(k-1)*m;
Ld[i] := L(n,m)
     end do
end do:
for i from 1 to N do
  for j from 1 to N do
    LL[i, j] := eval(Ld[i], t = T[j])
  end do
end do:

 
pn:=proc(i,n)
  if n=1 then
    return int(Ld[i],t)  
  else
    return int(pn(i,n-1),t)  
  end if
end proc:
##
RHS:= t->f(t)-a(t)*beta-b(t)*(t*beta+alpha):
f:=t->0:
a:=t->1/t:
b:=t->1:
alpha:=1:
beta:=0:
 

R:=Vector(N):
TMP:=Matrix(N,N):
C:=Matrix(N,N):
for i from 1 to N do
  R[i] := evalf(RHS(T[i]));
  tmp := a(t)*pn(i,1)+b(t)*pn(i,2);
  for j from 1 to N do
    TMP[i,j]:=eval(tmp, t = T[j])
  end do
end do:
use LinearAlgebra in
  C := Transpose(LinearSolve(Transpose(LL+TMP), R))
end use:
sol := beta*t+alpha+add(C[m0]*pn(m0,2), m0 = 1 .. N):
y:=unapply(sol,t):
ode:=t*diff(Y(t),t,t)+diff(Y(t),t)+t*Y(t)=0;  
 mExact:=unapply(rhs(dsolve({ode,Y(0)=1,D(Y)(0)=0})),t);
#
# Solve system using Maple numeric methods
# with all options on defaults
#
  mNumer:=eval( Y(t),
                dsolve
                ( {ode,Y(0)=1,D(Y)(0)=0},
                  numeric,
                  output=listprocedure
                )
              );

t*(diff(diff(Y(t), t), t))+diff(Y(t), t)+t*Y(t) = 0

 

proc (t) options operator, arrow; BesselJ(0, t) end proc

 

proc (t) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _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](t) else _xout := evalf(t) 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 .. 26, [( 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..63, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.10095317511683091e-1, (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..2, {(1) = 1.0, (2) = .0}, 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..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8]), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..2, {(1) = 1.0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = -.5}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..2, {(1, 1) = .0, (1, 2) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, t, Y, YP) option `[Y[1] = Y(t), Y[2] = diff(Y(t),t)]`; if t = 0 then if abs(Y[2]) <= 0. then YP[1] := 0; YP[2] := -(1/2)*Y[1] else error "system with provided initial conditions is singular" end if else YP[1] := Y[2]; YP[2] := (-t*Y[1]-Y[2])/t end if; 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, t, Y, YP) option `[Y[1] = Y(t), Y[2] = diff(Y(t),t)]`; if t = 0 then if abs(Y[2]) <= 0. then YP[1] := 0; YP[2] := -(1/2)*Y[1] else error "system with provided initial conditions is singular" end if else YP[1] := Y[2]; YP[2] := (-t*Y[1]-Y[2])/t end if; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0)  ] ))  ] ); _y0 := Array(0..2, {(1) = 0., (2) = 1.}); _vmap := array( 1 .. 2, [( 1 ) = (1), ( 2 ) = (2)  ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then error "initial conditions cannot be changed for systems with removable singularities"; _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); _i := false; if _par <> [] then _i := `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then _i := `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) or _i end if; if _i then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch:  end try else try procname(procname("right")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch:  end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch:  end try end if end if end if; return  elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error  end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; if type(_EnvDSNumericSaveDigits, 'posint') then _dat[4][26] := _EnvDSNumericSaveDigits else _dat[4][26] := Digits end if; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(1..3, {(1) = 18446744074581343590, (2) = 18446744074581343854, (3) = 18446744074581344030}), (3) = [t, Y(t), diff(Y(t), t)], (4) = []}); _solnproc := _dat[1]; _pars := map(rhs, _dat[4]); if not type(_xout, 'numeric') then if member(t, ["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(t, '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(t, ["last", 'last', "initial", 'initial', NULL]) then _res := _solnproc(convert(t, 'string')); if type(_res, 'list') then return _res[2] else return NULL end if elif member(t, ["parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(t, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[2], seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(t), 'string') = rhs(t); if lhs(_xout) = "initial" then if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else _res := _solnproc("initial" = ["single", 2, rhs(_xout)]) end if elif not type(rhs(_xout), 'list') then error "initial and/or parameter values must be specified in a list" elif lhs(_xout) = "initial_and_parameters" and nops(rhs(_xout)) = nops(_pars)+1 then _res := _solnproc(lhs(_xout) = ["single", 2, op(rhs(_xout))]) else _res := _solnproc(_xout) end if; if lhs(_xout) = "initial" then return _res[2] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[2], 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(t), 'string') = rhs(t)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _dat[3] end if; if procname <> unknown then return ('procname')(t) else _ndsol := 1; _ndsol := `tools/gensym`("Y(t)"); eval(FromInert(_Inert_FUNCTION(_Inert_NAME("assign"), _Inert_EXPSEQ(ToInert(_ndsol), _Inert_VERBATIM(pointto(_dat[2][2])))))); return FromInert(_Inert_FUNCTION(ToInert(_ndsol), _Inert_EXPSEQ(ToInert(t)))) end if end if; try _res := _solnproc(_xout); _res[2] catch: error  end try end proc

(1)

#
# Print all results to screen.
#
# Exact Maple solution
# Numerical Maple solution
# OP's Numerical solution
#
# Errors from Maple's numerical solution are slightly
# less than those obtained from OP's numerical method
#

  printf( "\n\n\n      Time    Maple Exact      "
          "  Maple Numeric       Maple Error    "
          "  OP Numeric          OP ERROR \n\n");
  seq( printf( " %8a %18.15f %18.15f %18.15f %18.15f %18.15f \n",
                 k,
                 evalf[12](mExact(k)),
                 mNumer(k),
                 mNumer(k)-mExact(k),
                 y(k),
                 evalf[12](mExact(k))-y(k) ),
       k=0..1, 1/10
     );




      Time    Maple Exact        Maple Numeric       Maple Error      OP Numeric          OP ERROR

        0  1.000000000000000  1.000000000000000  0.000000000000000  1.000000000000000  0.000000000000000
     1/10  0.997501562066000  0.997501562066904 -0.000000000033096  0.997501576700000 -0.000000014599130
      1/5  0.990024972240000  0.990024970551038 -0.000000001648962  0.990024978300000 -0.000000006101050
     3/10  0.977626246538000  0.977626243360794 -0.000000003139206  0.977626250400000 -0.000000003864050
      2/5  0.960398226660000  0.960398229592673  0.000000002892673  0.960398229600000 -0.000000002925130
      1/2  0.938469807241000  0.938469808879848  0.000000001679848  0.938469801300000  0.000000005880000
      3/5  0.912004863497000  0.912004845584140 -0.000000017915860  0.912004921400000 -0.000000057903760
     7/10  0.881200888607000  0.881200895274771  0.000000006674771  0.881200993000000 -0.000000104400960
      4/5  0.846287352750000  0.846287315917441 -0.000000036882559  0.846287497900000 -0.000000145100960
     9/10  0.807523798123000  0.807523793955245 -0.000000004144755  0.807523978900000 -0.000000180833760
        1  0.765197686558000  0.765197623203748 -0.000000063396252  0.765197891200000 -0.000000204600000

 

#
# Plot errors produced by Maple's numeric method
# and those produced by OP's numeric method
#
   OPErrors:=[seq([k,mExact(k)-y(k)], k=0..1, 1/10)]:
   MapleErrors:= [seq([k,mExact(k)-mNumer(k)], k=0..1, 1/10)]:
   plot( [OPErrors, MapleErrors], color =[red, blue], legend=["OP", "Maple"], axes=boxed);

 


Exact:=[seq([k,mExact(k)], k=0..1, 1/10)]:
Legendre:=[seq([k,y(k)], k=0..1, 1/10)]:
Error:=[seq([k,mExact(k)-y(k)], k=0..1, 1/10)]:
plot([Exact, Legendre, Error], color=[green,blue,red], legend=["Exact", "Legendre", "Error"], axes=box)

 

 


 

Download ODEnumer.mw

Use the menu setting

Tools->Options->Precision

and uncheck the box labelled "Limit expression length to"

I recomment that you use this option with care.

I would display the output here, but this site won't render it properly.

In a Maple worksheet I get a reasonably compact display of 10X6 graphs which fits on a single screen
 

in the following

  1. I have removed a lot of "stuff" which you don't really need
  2. In order to produce numeric values ( essential for plotting), I have assumed that 'm' is the electronic mass and 'hbar' is Planck's constant divided by 2*PI
  3. The obtained solution is complex (ie it contains I=sqrt(-1)). In ordert o guarantee evaluation to numeric vaues, the attached shows plots of abs( f(x,y,3) ), Re( f(x, y, 3) ) and Im( f(x, y, 3) )

Schrdiff(o(t), t, t)dinger PDE on (x,y,t) with initial and boundary conditions. Zero potential: problem number 169

 

Here is the problem 169 specification and solution from the Nasser list:

  restart;
  interface(rtablesize=10):
  interface(version);
  Physics:-Version();

`Standard Worksheet Interface, Maple 2019.0, Windows 7, March 9 2019 Build ID 1384062`

 

"C:\Users\TomLeslie\maple\toolbox\2019\Physics Updates\lib\Physics Updates.maple", `2019, April 26, 7:48 hours, version in the MapleCloud: 354, version installed in this computer: 350.`

(1.1)

  pde:= I*diff(f(x,y,t),t)=-hbar^2*(diff(f(x,y,t),x$2)+diff(f(x,y,t),y$2))/(2*m);
  ic:= f(x,y,0)=sqrt(2)*(sin(2*Pi*x)*sin(Pi*y)+sin(Pi*x)*sin(2*Pi*y));
  bc:= f(0,y,t)=0,
       f(1,y,t)=0,
       f(x,1,t)=0,
       f(x,0,t)=0;
  sol:= pdsolve([pde,ic,bc]);

I*(diff(f(x, y, t), t)) = -(1/2)*hbar^2*(diff(diff(f(x, y, t), x), x)+diff(diff(f(x, y, t), y), y))/m

 

f(x, y, 0) = 2^(1/2)*(sin(2*Pi*x)*sin(Pi*y)+sin(Pi*x)*sin(2*Pi*y))

 

f(0, y, t) = 0, f(1, y, t) = 0, f(x, 1, t) = 0, f(x, 0, t) = 0

 

f(x, y, t) = 2^(1/2)*sin(Pi*x)*exp(-((5/2)*I)*hbar^2*t*Pi^2/m)*(2*sin(Pi*y)*cos(Pi*x)+sin(2*Pi*y))

(1.2)

#
# Note: above solution is complex although
# the imaginary part is very small (see
# third plot)
#
# Subsequent plots give abs(), Re() and Im()
#
  plot3d
  ( eval
    ( abs(rhs(sol)),
      [ m = 9.10938356*10^(-30),
        hbar = 1.054571800*10^(-34),
        t = 3
      ]
    ),
    x = -1..1,
    y = -1..1,
    title= "abs( f( x, y, 3) )",
    titlefont= [ times, bold, 20]
  );
 plot3d
  ( eval
    ( Re(rhs(sol)),
      [ m = 9.10938356*10^(-30),
        hbar = 1.054571800*10^(-34),
        t = 3
      ]
    ),
    x = -1..1,
    y = -1..1,
    title= "Re( f( x, y, 3) )",
    titlefont= [ times, bold, 20]
  );
 plot3d
  ( eval
    ( Im(rhs(sol)),
      [ m = 9.10938356*10^(-30),
        hbar = 1.054571800*10^(-34),
        t = 3
      ]
    ),
    x = -1..1,
    y = -1..1,
    title= "Im( f( x, y, 3) )",
    titlefont= [ times, bold, 20]
  );

 

 

 

 

Download schPDE.mw

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