reduce-historical/r38/packages/odesolve/zimopbas.rlg

REDUCE Development Version, Wed Sep 13 20:40:41 2000 ...







ODESolve 1.065

% -*- REDUCE -*-
% The Postel/Zimmermann (11/4/96) ODE test examples.
% Equation names from Postel/Zimmermann.
% This version uses Maple-style functional notation wherever possible.
% It outputs general solutions of linear ODEs in basis format.
% It also checks all solutions.

on odesolve_basis, odesolve_check;


on div, intstr;

  off allfac;

            % to look prettier

% 1  Single equations without initial conditions
% ==============================================

% 1.1 Linear equations
% ====================

operator y;



% (1) Linear Bernoulli 1
odesolve((x^4-x^3)*df(y(x),x) + 2*x^4*y(x) = x^3/3 + C, y(x), x);


       - 2*x
     e
{{--------------},
    2
   x  - 2*x + 1

   1     -2    1       1
  ---*c*x   + ---*x - ---
   2           6       4
 -------------------------}
        2
       x  - 2*x + 1


% (2) Linear Bernoulli 2
odesolve(-1/2*df(y(x),x) + y(x) = sin x, y(x), x);


   2*x   2            4
{{e   },---*cos(x) + ---*sin(x)}
         5            5


% (3) Linear change of variables (FJW: shifted Euler equation)
odesolve(df(y(x),x,2)*(a*x+b)^2 + 4df(y(x),x)*(a*x+b)*a + 2y(x)*a^2 = 0,
   y(x), x);


             2
            a                a
{{----------------------,---------}}
    2  2              2   a*x + b
   a *x  + 2*a*b*x + b


% (4) Adjoint
odesolve((x^2-x)*df(y(x),x,2) + (2x^2+4x-3)*df(y(x),x) + 8x*y(x) = 1,
   y(x), x);


               2                                  2
{df(y(x),x,2)*x  - df(y(x),x,2)*x + 2*df(y(x),x)*x  + 4*df(y(x),x)*x

  - 3*df(y(x),x) + 8*y(x)*x - 1=0}


% (5) Polynomial solutions
% (FJW: Currently very slow, and fails anyway!)
% odesolve((x^2-x)*df(y(x),x,2) + (1-2x^2)*df(y(x),x) + (4x-2)*y(x) = 0,
%    y(x), x);

% (6) Dependent variable missing
odesolve(df(y(x),x,2) + 2x*df(y(x),x) = 2x, y(x), x);


   1
{{---*sqrt(pi)*erf(x),1},x}
   2


% (7) Liouvillian solutions
% (FJW: INTEGRATION IMPOSSIBLY SLOW WITHOUT EITHER ALGINT OR NOINT OPTION)
begin scalar !*allfac;  !*allfac := t;  return
   odesolve((x^3/2-x^2)*df(y(x),x,2) + (2x^2-3x+1)*df(y(x),x) + (x-1)*y(x) = 0,
      y(x), x, algint);
end;


               -1                         1/x
    - 1/2   - x           - 1/2          e   *sqrt(x - 2)
{{x      *e      *(x - 2)      *int(--------------------------,x),
                                              2
                                     sqrt(x)*x  - 2*sqrt(x)*x

               -1
    - 1/2   - x           - 1/2
  x      *e      *(x - 2)      }}

% NB: DO NOT RE-EVALUATE RESULT WITHOUT TURNING ON ALGINT OR NOINT SWITCH

% (8) Reduction of order
% (FJW: Attempting to make explicit currently too slow.)
odesolve(df(y(x),x,2) - 2x*df(y(x),x) + 2y(x) = 3, y(x), x);


***** Cannot convert nonlinear combination solution to basis!

{arbconst(2) + sqrt(pi)*arbconst(1)

                               erf(i*x)
 *int(-----------------------------------------------------------,x)*i
                                             2
                                            x
       sqrt(pi)*arbconst(1)*erf(i*x)*i*x + e  *arbconst(1) - 2*x

                                       1
  - 2*int(-----------------------------------------------------------,x)
                                                 2
                                                x
           sqrt(pi)*arbconst(1)*erf(i*x)*i*x + e  *arbconst(1) - 2*x

                3
  - log(y(x) - ---)=0}
                2


% (9) Integrating factors
% (FJW: Currently very slow, and fails anyway!)
% odesolve(sqrt(x)*df(y(x),x,2) + 2x*df(y(x),x) + 3y(x) = 0, y(x), x);

% (10) Radical solution (FJW: omitted for now)

% (11) Undetermined coefficients
odesolve(df(y(x),x,2) - 2/x^2*y(x) = 7x^4 + 3*x^3, y(x), x);


   -1  2   1   6    1   5
{{x  ,x },---*x  + ---*x }
           4        6


% (12) Variation of parameters
odesolve(df(y(x),x,2) + y(x) = csc(x), y(x), x);


{{cos(x),sin(x)}, - cos(x)*x + log(sin(x))*sin(x)}


% (13) Linear constant coefficients
<< factor exp(x);  write
odesolve(df(y(x),x,7) - 14df(y(x),x,6) + 80df(y(x),x,5) - 242df(y(x),x,4)
   + 419df(y(x),x,3) - 416df(y(x),x,2) + 220df(y(x),x) - 48y(x) = 0, y(x), x);
remfac exp(x) >>;


   3*x
{{e   ,

   4*x
  e   ,

   2*x
  e   *x,

   2*x
  e   ,

   x  2
  e *x ,

   x
  e *x,

   x
  e }}


% (14) Euler
odesolve(df(y(x),x,4) - 4/x^2*df(y(x),x,2) + 8/x^3*df(y(x),x) - 8/x^4*y(x) = 0,
   y(x), x);


   -1    2  4
{{x  ,x,x ,x }}


% (15) Exact n-th order
odesolve((1+x+x^2)*df(y(x),x,3) + (3+6x)*df(y(x),x,2) + 6df(y(x),x) = 6x,
   y(x), x);


      1   2
     ---*x
      2
{{------------,
    2
   x  + x + 1

       x
  ------------,
    2
   x  + x + 1

       1
  ------------},
    2
   x  + x + 1

     1   4
    ---*x
     4
 ------------}
   2
  x  + x + 1



% 1.2 Nonlinear equations
% =======================

% (16) Integrating factors 1
odesolve(df(y(x),x) = y(x)/(y(x)*log y(x) + x), y(x), x);


                       1           2
{x=arbconst(4)*y(x) + ---*log(y(x)) *y(x)}
                       2


% (17) Integrating factors 2
odesolve(2y(x)*df(y(x),x)^2 - 2x*df(y(x),x) - y(x) = 0, y(x), x);


                       4                 2      - 1/3             - 2/3  1/3
{{y(x)=2*(4*arbparam(1)  - 12*arbparam(1)  + 9)      *arbparam(1)      *2

  *arbconst(5)*arbparam(1),

                    4                 2      - 1/3             - 2/3  1/3
  x=2*(4*arbparam(1)  - 12*arbparam(1)  + 9)      *arbparam(1)      *2

                          2                 4                 2      - 1/3
  *arbconst(5)*arbparam(1)  - (4*arbparam(1)  - 12*arbparam(1)  + 9)

               - 2/3  1/3
  *arbparam(1)      *2   *arbconst(5),

  arbparam(1)}}

% This parametric solution is correct, cf. Zwillinger (1989) p.168 (41.10)
% (except that first edition is missing the constant C)!

% (18) Bernoulli 1
odesolve(df(y(x),x) + y(x) = y(x)^3*sin x, y(x), x, explicit);


{y(x)

      2*x                                    - 1/2
 =(5*e   *arbconst(6) + 2*cos(x) + 4*sin(x))      *sqrt(5)*plus_or_minus(tag_1)}

expand_plus_or_minus ws;


          2*x                                    - 1/2
{y(x)=(5*e   *arbconst(6) + 2*cos(x) + 4*sin(x))      *sqrt(5),

             2*x                                    - 1/2
 y(x)= - (5*e   *arbconst(6) + 2*cos(x) + 4*sin(x))      *sqrt(5)}


% (19) Bernoulli 2
operator P, Q;


begin scalar soln, !*exp, !*allfac;  % for a neat solution
   on allfac;
   soln := odesolve(df(y(x),x) + P(x)*y(x) = Q(x)*y(x)^n, y(x), x);
   off allfac;  return soln
end;


      - n          (n - 1)*int(p(x),x)
{y(x)    *y(x)= - e

                 - int(p(x),x)*n + int(p(x),x)
 *((n - 1)*int(e                              *q(x),x) - arbconst(7))}

odesolve(df(y(x),x) + P(x)*y(x) = Q(x)*y(x)^(2/3), y(x), x);


     1/3   - 1/3*int(p(x),x)
{y(x)   =e                  *arbconst(8)

     1    - 1/3*int(p(x),x)      1/3*int(p(x),x)
  + ---*e                  *int(e               *q(x),x)}
     3


% (20) Clairaut 1
odesolve((x^2-1)*df(y(x),x)^2 - 2x*y(x)*df(y(x),x) + y(x)^2 - 1 = 0,
   y(x), x, explicit);


                                      2
{y(x)=arbconst(9)*x + sqrt(arbconst(9)  + 1),

                                      2
 y(x)=arbconst(9)*x - sqrt(arbconst(9)  + 1),

               2
 y(x)=sqrt( - x  + 1),

                  2
 y(x)= - sqrt( - x  + 1)}


% (21) Clairaut 2
operator f, g;


odesolve(f(x*df(y(x),x)-y(x)) = g(df(y(x),x)), y(x), x);


{f(arbconst(10)*x - y(x)) - g(arbconst(10))=0}


% (22) Equations of the form  y' = f(x,y)
odesolve(df(y(x),x) = (3x^2-y(x)^2-7)/(exp(y(x))+2x*y(x)+1), y(x), x);


                 y(x)       2             3
{arbconst(11) + e     + y(x) *x + y(x) - x  + 7*x=0}


% (23) Homogeneous
odesolve(df(y(x),x) = (2x^3*y(x)-y(x)^4)/(x^4-2x*y(x)^3), y(x), x);


                           3    3
{arbconst(12)*y(x)*x + y(x)  + x =0}


% (24) Factoring the equation
odesolve(df(y(x),x)*(df(y(x),x)+y(x)) = x*(x+y(x)), y(x), x);


        - x
{y(x)=e    *arbconst(13) - x + 1,

                      1   2
 y(x)=arbconst(14) + ---*x }
                      2


% (25) Interchange variables
% (NB: Soln in Zwillinger (1989) wrong, as is last eqn in Table 68!)
odesolve(df(y(x),x) = x/(x^2*y(x)^2+y(x)^5), y(x), x);


             3
  2  2/3*y(x)                     3    3
{x =e         *arbconst(15) - y(x)  - ---}
                                       2


% (26) Lagrange 1
odesolve(y(x) = 2x*df(y(x),x) - a*df(y(x),x)^3, y(x), x);


                                 -1    1             3
{{y(x)=2*arbconst(16)*arbparam(2)   + ---*arbparam(2) *a,
                                       2

                            -2    3             2
  x=arbconst(16)*arbparam(2)   + ---*arbparam(2) *a,
                                  4

  arbparam(2)}}

odesolve(y(x) = 2x*df(y(x),x) - a*df(y(x),x)^3, y(x), x, implicit);


                3  2                   2    2                        2
{64*arbconst(17) *a  + 128*arbconst(17) *a*x  - 144*arbconst(17)*y(x) *a*x

                     4          4            2  3
  + 64*arbconst(17)*x  + 27*y(x) *a - 16*y(x) *x =0}

% root_of quartic is VERY slow if explicit option used!

% (27) Lagrange 2
odesolve(y(x) = 2x*df(y(x),x) - df(y(x),x)^2, y(x), x);


                                 -1    1             2
{{y(x)=2*arbconst(18)*arbparam(3)   + ---*arbparam(3) ,
                                       3

                            -2    2
  x=arbconst(18)*arbparam(3)   + ---*arbparam(3),
                                  3

  arbparam(3)}}

odesolve(y(x) = 2x*df(y(x),x) - df(y(x),x)^2, y(x), x, implicit);


                  2                                             3         3
{ - 9*arbconst(19)  + 18*arbconst(19)*y(x)*x - 12*arbconst(19)*x  - 4*y(x)

          2  2
  + 3*y(x) *x =0}


% (28) Riccati 1
odesolve(df(y(x),x) = exp(x)*y(x)^2 - y(x) + exp(-x), y(x), x);


         - x                         - x
       e    *arbconst(20)*sin(x) - e    *cos(x)
{y(x)=------------------------------------------}
             arbconst(20)*cos(x) + sin(x)


% (29) Riccati 2
<< factor x;  write
odesolve(df(y(x),x) = y(x)^2 - x*y(x) + 1, y(x), x);
remfac x >>;


                                  2
                             1/2*x
                          2*e      *arbconst(21)
{y(x)=x + ------------------------------------------------------}
                                               - 1/2
           sqrt(pi)*sqrt(2)*arbconst(21)*erf(2      *x*i)*i - 2


% (30) Separable
odesolve(df(y(x),x) = (9x^8+1)/(y(x)^2+1), y(x), x);


                      3               9
{3*arbconst(22) + y(x)  + 3*y(x) - 3*x  - 3*x=0}


% (31) Solvable for x
odesolve(y(x) = 2x*df(y(x),x) + y(x)*df(y(x),x)^2, y(x), x);


                                    -1
{{y(x)= - 2*arbconst(23)*arbparam(4)  ,

                               -2
  x= - arbconst(23)*arbparam(4)   + arbconst(23),

  arbparam(4)}}

odesolve(y(x) = 2x*df(y(x),x) + y(x)*df(y(x),x)^2, y(x), x, implicit);


                  2                          2
{ - 4*arbconst(24)  + 4*arbconst(24)*x + y(x) =0}


% (32) Solvable for y
begin scalar !*allfac;  !*allfac := t;  return
   odesolve(x = y(x)*df(y(x),x) - x*df(y(x),x)^2, y(x), x)
end;


                          2
         - 1/2*arbparam(5)                           2
{{y(x)=e                   *arbconst(25)*(arbparam(5)  + 1),

                       2
      - 1/2*arbparam(5)
  x=e                   *arbconst(25)*arbparam(5),

  arbparam(5)}}


% (33) Autonomous 1
odesolve(df(y(x),x,2)-df(y(x),x) = 2y(x)*df(y(x),x), y(x), x, explicit);


{y(x)

      1                    arbconst(27)*arbconst(26) - arbconst(26)*x      1
 = - ---*arbconst(26)*tan(--------------------------------------------) - ---,
      2                                        2                           2

 y(x)=arbconst(28)}


% (34) Autonomous 2  (FJW: Slow without either algint or noint option.)
odesolve(df(y(x),x,2)/y(x) - df(y(x),x)^2/y(x)^2 - 1 + 1/y(x)^3 = 0,
   y(x), x, algint);


{arbconst(31)*plus_or_minus(tag_4) + sqrt(3)

                                     3                   3      - 1/2
 *int(sqrt(y(x))*(3*arbconst(30)*y(x)  + 6*log(y(x))*y(x)  + 2)      ,y(x))

  - plus_or_minus(tag_4)*x=0}


% (35) Differentiation method
odesolve(2y(x)*df(y(x),x,2) - df(y(x),x)^2 =
   1/3(df(y(x),x) - x*df(y(x),x,2))^2, y(x), x, explicit);


                  2  2                                                         2
{y(x)=arbconst(33) *x  + 2*sqrt(3)*arbconst(33)*arbconst(32)*x + 4*arbconst(32)

 ,

                  2  2                                                         2
 y(x)=arbconst(34) *x  - 2*sqrt(3)*arbconst(34)*arbconst(32)*x + 4*arbconst(32)

 ,

 y(x)=arbconst(35)}


% (36) Equidimensional in x
odesolve(x*df(y(x),x,2) = 2y(x)*df(y(x),x), y(x), x, explicit);


          1                    arbconst(37)*arbconst(36) - arbconst(36)*log(x)
{y(x)= - ---*arbconst(36)*tan(-------------------------------------------------)
          2                                           2

     1
  - ---,
     2

 y(x)=arbconst(38)}


% (37) Equidimensional in y
odesolve((1-x)*(y(x)*df(y(x),x,2)-df(y(x),x)^2) + x^2*y(x)^2 = 0, y(x), x);


                                             3        2
        arbconst(40) + arbconst(39)*x + 1/6*x  + 1/2*x  - x        x
       e                                                   *(x - 1)
{y(x)=---------------------------------------------------------------}
                                   x - 1


% (38) Exact second order
odesolve(x*y(x)*df(y(x),x,2) + x*df(y(x),x)^2 + y(x)*df(y(x),x) = 0,
   y(x), x, explicit);


{y(x)=sqrt( - arbconst(42) + log(x))*sqrt(arbconst(41))*sqrt(2),

 y(x)= - sqrt( - arbconst(42) + log(x))*sqrt(arbconst(41))*sqrt(2),

 y(x)=arbconst(43)}


% (39) Factoring differential operator
odesolve(df(y(x),x,2)^2 - 2df(y(x),x)*df(y(x),x,2) + 2y(x)*df(y(x),x) -
   y(x)^2 = 0, y(x), x);


       x                 x
{y(x)=e *arbconst(45) + e *arbconst(44)*x,

       x                  - x
 y(x)=e *arbconst(47) + e    *arbconst(46)}


% (40) Scale invariant (fails with algint option)
odesolve(x^2*df(y(x),x,2) + 3x*df(y(x),x) = 1/(y(x)^3*x^4), y(x), x);


{2*arbconst(49)*plus_or_minus(tag_7) + log(

                        2      - 1/2                               2      - 1/2
     - 2*(4*arbconst(48)  + 1)      *arbconst(48) + (4*arbconst(48)  + 1)

                                2  2       4  4
    *sqrt( - 4*arbconst(48)*y(x) *x  + y(x) *x  - 1)

                      2      - 1/2     2  2
     + (4*arbconst(48)  + 1)      *y(x) *x ) - 2*log(x)*plus_or_minus(tag_7)=0}


% Revised scale-invariant example (hangs with algint option):
ode := x^2*df(y(x),x,2) + 3x*df(y(x),x) + 2*y(x) = 1/(y(x)^3*x^4);


                     2                               -3  -4
ode := df(y(x),x,2)*x  + 3*df(y(x),x)*x + 2*y(x)=y(x)  *x

% Choose full (explicit and expanded) solution:
odesolve(ode, y(x), x, full);


          1
{y(x)= - ---*sqrt(15*arbconst(50)
          2

                            2                                         - 1/2  -1
     - sqrt(225*arbconst(50)  - 64)*sin(2*(arbconst(51) - log(x))))*2      *x  ,

          1
 y(x)= - ---*sqrt(15*arbconst(50)
          2

                            2                                         - 1/2  -1
     + sqrt(225*arbconst(50)  - 64)*sin(2*(arbconst(51) - log(x))))*2      *x  ,

       1
 y(x)=---*sqrt(15*arbconst(50)
       2

                                       2
                - sqrt(225*arbconst(50)  - 64)*sin(2*(arbconst(51) - log(x))))

    - 1/2  -1
 *2      *x  ,

       1
 y(x)=---*sqrt(15*arbconst(50)
       2

                                       2
                + sqrt(225*arbconst(50)  - 64)*sin(2*(arbconst(51) - log(x))))

    - 1/2  -1
 *2      *x  }
           % or "explicit, expand"
% Check it -- each solution should simplify to 0:
foreach soln in ws collect
   trigsimp sub(soln, num(lhs ode - rhs ode));


{0,0,0,0}


% (41) Autonomous, 3rd order
odesolve((df(y(x),x)^2+1)*df(y(x),x,3) - 3df(y(x),x)*df(y(x),x,2)^2 = 0,
   y(x), x);


                                      2             2
{y(x)=arbconst(55) + sqrt(arbconst(53) *arbconst(52)

                                  2                                  2  2
     - 2*arbconst(53)*arbconst(52) *x + 2*arbconst(53) + arbconst(52) *x  - 2*x)

              -1
 *arbconst(52)  *i,

 y(x)=arbconst(56) + i*x,

 y(x)=arbconst(57) - i*x,

 y(x)=arbconst(58) + arbconst(54)*x}

% odesolve((df(y(x),x)^2+1)*df(y(x),x,3) - 3df(y(x),x)*df(y(x),x,2)^2 = 0,
%    y(x), x, implicit);
% Implicit form is currently too messy!

% (42) Autonomous, 4th order
odesolve(3*df(y(x),x,2)*df(y(x),x,4) - 5df(y(x),x,3)^2 = 0, y(x), x);


{y(x)=arbconst(63)*x + arbconst(62)

                                                 -3
  - 3*sqrt(arbconst(60) - x)*sqrt(6)*arbconst(59)  ,

 y(x)=arbconst(65)*x + arbconst(64)

                                                 -3
  + 3*sqrt(arbconst(60) - x)*sqrt(6)*arbconst(59)  ,

                                       1                2
 y(x)=arbconst(67)*x + arbconst(66) + ---*arbconst(61)*x }
                                       2



% 1.3 Special equations
% =====================

% (43) Delay
odesolve(df(y(x),x) + a*y(x-1) = 0, y(x), x);


***** Arguments of y differ -- solving delay equations is not implemented.


% (44) Functions with several parameters
odesolve(df(y(x,a),x) = a*y(x,a), y(x,a), x);


   a*x
{{e   }}



% 2  Single equations with initial conditions
% ===========================================

% (45) Exact 4th order
odesolve(df(y(x),x,4) = sin x, y(x), x,
   {x=0, y(x)=0, df(y(x),x)=0, df(y(x),x,2)=0, df(y(x),x,3)=0});


                1   3
{y(x)=sin(x) + ---*x  - x}
                6


% (46) Linear polynomial coefficients -- Bessel J0
odesolve(x*df(y(x),x,2) + df(y(x),x) + 2x*y(x) = 0, y(x), x,
   {x=0, y(x)=1, df(y(x),x)=0});


{y(x)=besselj(0,sqrt(2)*x)}


% (47) Second-degree separable
soln :=
   odesolve(x*df(y(x),x)^2 - y(x)^2 + 1 = 0, y(x)=1, x=0, explicit);


               1   2*sqrt(x)*plus_or_minus(tag_9)
soln := {y(x)=---*e
               2

             1    - 2*sqrt(x)*plus_or_minus(tag_9)
          + ---*e                                 }
             2

% Alternatively ...
soln where e^~x => cosh x + sinh x;


{y(x)=cosh(2*sqrt(x)*plus_or_minus(tag_9))}

% but this works ONLY with `on div, intstr; off allfac;'
% A better alternative is ...
trigsimp(soln, hyp, combine);


{y(x)=cosh(2*sqrt(x)*plus_or_minus(tag_9))}

expand_plus_or_minus ws;


{y(x)=cosh(2*sqrt(x))}


% (48) Autonomous
odesolve(df(y(x),x,2) + y(x)*df(y(x),x)^3 = 0, y(x), x,
   {x=0, y(x)=0, df(y(x),x)=2});


     3
{y(x)  + 3*y(x) - 6*x=0}

%% Only one explicit solution satisfies the conditions:
begin scalar !*trode, !*fullroots;  !*fullroots := t;  return
   odesolve(df(y(x),x,2) + y(x)*df(y(x),x)^3 = 0, y(x), x,
      {x=0, y(x)=0, df(y(x),x)=2}, explicit);
end;


               2            1/3            2             - 1/3
{y(x)=(sqrt(9*x  + 1) + 3*x)    - (sqrt(9*x  + 1) + 3*x)      }



% 3  Systems of equations
% =======================

% (49) Integrable combinations
operator x, z;


odesolve({df(x(t),t) = -3y(t)*z(t), df(y(t),t) = 3x(t)*z(t),
   df(z(t),t) = -x(t)*y(t)}, {x(t),y(t),z(t)}, t);


odesolve-system({df(x(t),t) + 3*y(t)*z(t),

                 df(y(t),t) - 3*x(t)*z(t),

                 df(z(t),t) + x(t)*y(t)},{x(t),y(t),z(t)},t)


% (50) Matrix Riccati
operator a, b;


odesolve({df(x(t),t) = a(t)*(y(t)^2-x(t)^2) + 2b(t)*x(t)*y(t) + 2c*x(t),
   df(y(t),t) = b(t)*(y(t)^2-x(t)^2) - 2a(t)*x(t)*y(t) + 2c*y(t)},
   {x(t),y(t)}, t);


                          2            2
odesolve-system({a(t)*x(t)  - a(t)*y(t)  - 2*b(t)*x(t)*y(t) + df(x(t),t)

                  - 2*c*x(t),

                                             2            2
                 2*a(t)*x(t)*y(t) + b(t)*x(t)  - b(t)*y(t)  + df(y(t),t)

                  - 2*c*y(t)},{x(t),y(t)},t)


% (51) Triangular
odesolve({df(x(t),t) = x(t)*(1 + cos(t)/(2+sin(t))),
   df(y(t),t) = x(t) - y(t)}, {x(t),y(t)}, t);


odesolve-system({( - cos(t)*x(t) + df(x(t),t)*sin(t) + 2*df(x(t),t)

                     - sin(t)*x(t) - 2*x(t))/(sin(t) + 2),

                 df(y(t),t) - x(t) + y(t)},{x(t),y(t)},t)


% (52) Vector
odesolve({df(x(t),t) = 9x(t) + 2y(t), df(y(t),t) = x(t) + 8y(t)},
   {x(t),y(t)}, t);


odesolve-system({df(x(t),t) - 9*x(t) - 2*y(t),

                 df(y(t),t) - x(t) - 8*y(t)},{x(t),y(t)},t)


% (53) Higher order
odesolve({df(x(t),t) - x(t) + 2y(t) = 0,
   df(x(t),t,2) - 2df(y(t),t) = 2t - cos(2t)}, {x(t),y(t)}, t);


odesolve-system({df(x(t),t) - x(t) + 2*y(t),

                 cos(2*t) + df(x(t),t,2) - 2*df(y(t),t) - 2*t},{x(t),y(t)},t)


% (54) Inhomogeneous system
equ := {df(x(t),t) = -1/(t*(t^2+1))*x(t) + 1/(t^2*(t^2+1))*y(t) + 1/t,
        df(y(t),t) = -t^2/(t^2+1)*x(t) + (2t^2+1)/(t*(t^2+1))*y(t) + 1};


                             -1         -2        -1
                     - x(t)*t   + y(t)*t   + t + t
equ := {df(x(t),t)=----------------------------------,
                                  2
                                 t  + 1

                             2                    -1    2
                     - x(t)*t  + 2*y(t)*t + y(t)*t   + t  + 1
        df(y(t),t)=-------------------------------------------}
                                      2
                                     t  + 1

odesolve(equ, {x(t),y(t)}, t);


                              2                     -1         -1         -2
                  df(x(t),t)*t  + df(x(t),t) - t + t  *x(t) - t   - y(t)*t
odesolve-system({------------------------------------------------------------,
                                             2
                                            t  + 1

                              2                 2         2
                 (df(y(t),t)*t  + df(y(t),t) + t *x(t) - t  - 2*t*y(t)

                             -1        2
                     - y(t)*t   - 1)/(t  + 1)},{x(t),y(t)},t)


end;



Time for test: 23953 ms, plus GC time: 1807 ms