close all
clear all
clc

% Symbolic Solution of ODEs
% Page numbers are from the instructor's notes.
% Ugur Akman 2004

% The function 'dsolve' computes symbolic solutions to ordinary differential
% equations. The equations are specified by symbolic expressions containing
% the letter 'D' to denote differentiation. The symbols D2, D3, ... DN,
% correspond to the second, third, ..., Nth derivative, respectively. Thus,
% D2y is the Symbolic Math Toolbox equivalent of d^2y/dt^2. The dependent
% variables are those preceded by D and the default independent variable
% is t. Note that names of symbolic variables should not contain D. The
% independent variable can be changed from t to some other symbolic variable
% by including that variable as the last input argument. Initial conditions
% can be specified by additional equations. If initial conditions are not
% specified, the solutions contain constants of integration, C1, C2, etc.
% The output from 'dsolve' parallels the output from 'solve'. That is, you
% can call dsolve with the number of output variables equal to the number of
% dependent variables or place the output in a structure whose fields
% contain the solutions of the differential equations.


% Example - p.19
ODE = 'Dy = -(x*(1+y^2)^(1/2)) / (y*(1+x^2)^(1/2))'
% or ODE = 'Dy*(y*(1+x^2)^(1/2)) + (x*(1+y^2)^(1/2)) = 0'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.20
ODE = '(x*y^2+y) - x*Dy = 0'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.21
ODE = 'Dy = (y-4*x)^2 '
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.23
ODE = '(x+y) + (3*x+3*y-4)*Dy = 0'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.24
ODE = '(x-y-1) + (4*y+x-1)*Dy = 0'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.25
ODE = '(1-x*y-(x^2)*(y^2)) + ((x^3)*y-x^2)*Dy = 0'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.27
ODE = 'x*Dy + y = 0'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.27
ODE = '(cos(x)+2*y)*Dy = -x^3 + y*sin(x)'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.28
ODE = 'Dy = -2*x*y + 4*x'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.31 without initial condition
ODE = 'DT + a*T = 10*a*(1+cos(b*t))'
Solution = dsolve(ODE, 't')
Simplified = simplify(Solution)
pretty(Solution)

% Example - p.31 with initial condition
ODE = 'DT + a*T = 10*a*(1+cos(b*t))'
Solution = dsolve(ODE, 'T(0)=T0', 't')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.32
ODE = 'x*Dy = y+x*y^3*(1-log(x))'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.34
ODE = '(x^2+x)*(Dy)^2 + (x^2+x-2*x*y-y)*Dy-y^2-x*y = 0'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.35
ODE = '(Dy)^2 = y - (2 + Dy)*x'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.41
ODE = 'D2y + x*Dy = x'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)

% Example - p.49
ODE = 'D3y - 3*D2y + 3*Dy = x'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)
%simple(Solution)

% Example - p.50
ODE = '3*D2y + 10*Dy - 8*y = 7*exp(-4*x)'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.52
ODE = 'D2y - 4*Dy + 4*y = 2*exp(2*x) + cos(x)'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.56
ODE = 'D2y + y = tan(x)'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.57
ODE = 'x*D2y + Dy = x+1'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.58
ODE = 'x*D2y - 2*x*Dy + 2*y = x^2+2'
Solution = dsolve(ODE, 'x')
Simplified = simplify(Solution)
pretty(Solution)


% Example - p.60
ODE1 = 'D2x - x - 2*y = t'
ODE2 = 'D2y - 2*y - 3*x = 1'
[Solution1 Solution2] = dsolve(ODE1, ODE2, 't')
Simplified1 = simplify(Solution1)
Simplified2 = simplify(Solution2)
pretty(Solution1)
pretty(Solution2)


% Example - p.62
ODE1 = 'Dx - y = t'
ODE2 = 'Dy - x = 1'
[Solution1 Solution2] = dsolve(ODE1, ODE2, 't')
Simplified1 = simplify(Solution1)
Simplified2 = simplify(Solution2)
pretty(Solution1)
pretty(Solution2)


% Example - p.62 with initial conditions & plot
ODE1 = 'Dx - y = t'
ODE2 = 'Dy - x = 1'
[Solution1 Solution2] = dsolve(ODE1, ODE2, 'x(0)=0.1', 'y(0)=0.1', 't')
t=0:0.01:1;
EvalSimp1=eval(Solution1);
EvalSimp2=eval(Solution2);
figure
plot(t,EvalSimp1, 'r')
hold on
plot(t,EvalSimp2, 'b')
ODE =

Dy = -(x*(1+y^2)^(1/2)) / (y*(1+x^2)^(1/2))

Warning: Explicit solution could not be found; implicit solution returned.
 
Solution =
 
(1+x^2)^(1/2)+(1+y^2)^(1/2)+C1 = 0
 
 
 
Simplified =
 
(1+x^2)^(1/2)+(1+y^2)^(1/2)+C1 = 0
 
 
 
                            2 1/2         2 1/2
                      (1 + x )    + (1 + y )    + C1 = 0

ODE =

(x*y^2+y) - x*Dy = 0

 
Solution =
 
-2*x/(x^2-2*C1)
 
 
 
Simplified =
 
-2*x/(x^2-2*C1)
 
 
 
                                        x
                                 -2 ---------
                                     2
                                    x  - 2 C1

ODE =

Dy = (y-4*x)^2 

 
Solution =
 
2*(-exp(x)^4*C1-2*x+2*x*exp(x)^4*C1-1)/(-1+exp(x)^4*C1)
 
 
 
Simplified =
 
2*(-exp(4*x)*C1-2*x+2*x*exp(4*x)*C1-1)/(-1+exp(4*x)*C1)
 
 
 
                            4                      4
                     -exp(x)  C1 - 2 x + 2 x exp(x)  C1 - 1
                   2 --------------------------------------
                                           4
                                -1 + exp(x)  C1

ODE =

(x+y) + (3*x+3*y-4)*Dy = 0

 
Solution =
 
exp(-lambertw(3/2*exp(x)*exp(-3)/exp(C1))+x-3-C1)-x+2
 
 
 
Simplified =
 
2/3*lambertw(3/2*exp(x-3-C1))-x+2
 
 
 
                              exp(x) exp(-3)
            exp(-lambertw(3/2 --------------) + x - 3 - C1) - x + 2
                                 exp(C1)

ODE =

(x-y-1) + (4*y+x-1)*Dy = 0

Warning: Explicit solution could not be found; implicit solution returned.
 
Solution =
 
-1/2*log(((x-1)^2+4*y^2)/(x-1)^2)+1/2*atan(-2*y/(x-1))-log(x-1)-C1 = 0
 
 
 
Simplified =
 
-1/2*log((x^2-2*x+1+4*y^2)/(x-1)^2)-1/2*atan(2*y/(x-1))-log(x-1)-C1 = 0
 
 
 
                      2      2
               (x - 1)  + 4 y                  y
     - 1/2 log(---------------) - 1/2 atan(2 -----) - log(x - 1) - C1 = 0
                         2                   x - 1
                  (x - 1)

ODE =

(1-x*y-(x^2)*(y^2)) + ((x^3)*y-x^2)*Dy = 0

Warning: Explicit solution could not be found; implicit solution returned.
 
Solution =
 
log(x)-C1-1/4*log(2*x^2*y^2-1)-1/2*2^(1/2)*atanh(x*y*2^(1/2)) = 0
 
 
 
Simplified =
 
log(x)-C1-1/4*log(2*x^2*y^2-1)-1/2*2^(1/2)*atanh(x*y*2^(1/2)) = 0
 
 
 
                                2  2             1/2            1/2
       log(x) - C1 - 1/4 log(2 x  y  - 1) - 1/2 2    atanh(x y 2   ) = 0

ODE =

x*Dy + y = 0

 
Solution =
 
C1/x
 
 
 
Simplified =
 
C1/x
 
 
 
                                      C1
                                     ----
                                      x

ODE =

(cos(x)+2*y)*Dy = -x^3 + y*sin(x)

 
Solution =
 
 -1/2*cos(x)-1/2*(cos(x)^2-x^4-4*C1)^(1/2)
 -1/2*cos(x)+1/2*(cos(x)^2-x^4-4*C1)^(1/2)
 
 
 
Simplified =
 
 -1/2*cos(x)-1/2*(cos(x)^2-x^4-4*C1)^(1/2)
 -1/2*cos(x)+1/2*(cos(x)^2-x^4-4*C1)^(1/2)
 
 
 
                 [                          2    4        1/2]
                 [- 1/2 cos(x) - 1/2 (cos(x)  - x  - 4 C1)   ]
                 [                                           ]
                 [                          2    4        1/2]
                 [- 1/2 cos(x) + 1/2 (cos(x)  - x  - 4 C1)   ]

ODE =

Dy = -2*x*y + 4*x

 
Solution =
 
2+exp(-x^2)*C1
 
 
 
Simplified =
 
2+exp(-x^2)*C1
 
 
 
                                          2
                                2 + exp(-x ) C1

ODE =

DT + a*T = 10*a*(1+cos(b*t))

 
Solution =
 
exp(-a*t)*C1+10*(a^2+b^2+a^2*cos(b*t)+a*b*sin(b*t))/(a^2+b^2)
 
 
 
Simplified =
 
(exp(-a*t)*C1*a^2+exp(-a*t)*C1*b^2+10*a^2+10*b^2+10*a^2*cos(b*t)+10*a*b*sin(b*t))/(a^2+b^2)
 
 
 
                               2    2    2
                              a  + b  + a  cos(b t) + a b sin(b t)
            exp(-a t) C1 + 10 ------------------------------------
                                             2    2
                                            a  + b

ODE =

DT + a*T = 10*a*(1+cos(b*t))

 
Solution =
 
exp(-a*t)*(-20*a^2-10*b^2+T0*a^2+T0*b^2)/(a^2+b^2)+10*(a^2+b^2+a^2*cos(b*t)+a*b*sin(b*t))/(a^2+b^2)
 
 
 
Simplified =
 
(-20*exp(-a*t)*a^2-10*exp(-a*t)*b^2+exp(-a*t)*T0*a^2+exp(-a*t)*T0*b^2+10*a^2+10*b^2+10*a^2*cos(b*t)+10*a*b*sin(b*t))/(a^2+b^2)
 
 
 
                  2       2       2       2
  exp(-a t) (-20 a  - 10 b  + T0 a  + T0 b )
  ------------------------------------------
                    2    2
                   a  + b

               2    2    2
              a  + b  + a  cos(b t) + a b sin(b t)
         + 10 ------------------------------------
                             2    2
                            a  + b

ODE =

x*Dy = y+x*y^3*(1-log(x))

 
Solution =
 
  3/(-8*x^3+6*x^3*log(x)+9*C1)^(1/2)*x
 -3/(-8*x^3+6*x^3*log(x)+9*C1)^(1/2)*x
 
 
 
Simplified =
 
  3/(-8*x^3+6*x^3*log(x)+9*C1)^(1/2)*x
 -3/(-8*x^3+6*x^3*log(x)+9*C1)^(1/2)*x
 
 
 
                     [                 x                ]
                     [3 ------------------------------- ]
                     [       3      3               1/2 ]
                     [  (-8 x  + 6 x  log(x) + 9 C1)    ]
                     [                                  ]
                     [                  x               ]
                     [-3 -------------------------------]
                     [        3      3               1/2]
                     [   (-8 x  + 6 x  log(x) + 9 C1)   ]

ODE =

(x^2+x)*(Dy)^2 + (x^2+x-2*x*y-y)*Dy-y^2-x*y = 0

Warning: Explicit solution could not be found.
 
Solution =
 
[ empty sym ]
 
 
 
Simplified =
 
[ empty sym ]
 
 
 
                                   array([])

ODE =

(Dy)^2 = y - (2 + Dy)*x

 
Solution =
 
(4-1/2*x+2*lambertw(1/4*C1*exp(1/4*x-1)))*x+(-1/2*x+2+2*lambertw(1/4*C1*exp(1/4*x-1)))^2
 
 
 
Simplified =
 
2*x-1/4*x^2+4+8*lambertw(1/4*C1*exp(1/4*x-1))+4*lambertw(1/4*C1*exp(1/4*x-1))^2
 
 
 
  (4 - 1/2 x + 2 lambertw(1/4 C1 exp(1/4 x - 1))) x

                                                            2
         + (- 1/2 x + 2 + 2 lambertw(1/4 C1 exp(1/4 x - 1)))

ODE =

D2y + x*Dy = x

 
Solution =
 
1/2*C1*pi^(1/2)*2^(1/2)*erf(1/2*x*2^(1/2))+x+C2
 
 
 
Simplified =
 
1/2*C1*pi^(1/2)*2^(1/2)*erf(1/2*x*2^(1/2))+x+C2
 
 
 
                           1/2  1/2            1/2
                  1/2 C1 pi    2    erf(1/2 x 2   ) + x + C2

ODE =

D3y - 3*D2y + 3*Dy = x

 
Solution =
 
-1/6*C2*3^(1/2)*exp(3/2*x)*cos(1/2*3^(1/2)*x)+1/2*exp(3/2*x)*sin(1/2*3^(1/2)*x)*C2+1/2*exp(3/2*x)*cos(1/2*3^(1/2)*x)*C1+1/6*C1*3^(1/2)*exp(3/2*x)*sin(1/2*3^(1/2)*x)+1/6*x^2+1/3*x+C3
 
 
 
Simplified =
 
-1/6*C2*3^(1/2)*exp(3/2*x)*cos(1/2*3^(1/2)*x)+1/2*exp(3/2*x)*sin(1/2*3^(1/2)*x)*C2+1/2*exp(3/2*x)*cos(1/2*3^(1/2)*x)*C1+1/6*C1*3^(1/2)*exp(3/2*x)*sin(1/2*3^(1/2)*x)+1/6*x^2+1/3*x+C3
 
 
 
            1/2                     1/2                              1/2
  - 1/6 C2 3    exp(3/2 x) cos(1/2 3    x) + 1/2 exp(3/2 x) sin(1/2 3    x) C2

                                   1/2
         + 1/2 exp(3/2 x) cos(1/2 3    x) C1

                   1/2                     1/2           2
         + 1/6 C1 3    exp(3/2 x) sin(1/2 3    x) + 1/6 x  + 1/3 x + C3

ODE =

3*D2y + 10*Dy - 8*y = 7*exp(-4*x)

 
Solution =
 
exp(-4*x)*C2+exp(2/3*x)*C1-1/2*x*exp(-4*x)
 
 
 
Simplified =
 
exp(-4*x)*C2+exp(2/3*x)*C1-1/2*x*exp(-4*x)
 
 
 
                exp(-4 x) C2 + exp(2/3 x) C1 - 1/2 x exp(-4 x)

ODE =

D2y - 4*Dy + 4*y = 2*exp(2*x) + cos(x)

 
Solution =
 
exp(2*x)*C2+exp(2*x)*x*C1+exp(2*x)*x^2+3/25*cos(x)-4/25*sin(x)
 
 
 
Simplified =
 
exp(2*x)*C2+exp(2*x)*x*C1+exp(2*x)*x^2+3/25*cos(x)-4/25*sin(x)
 
 
 
                                             2
     exp(2 x) C2 + exp(2 x) x C1 + exp(2 x) x  + 3/25 cos(x) - 4/25 sin(x)

ODE =

D2y + y = tan(x)

 
Solution =
 
sin(x)*C2+cos(x)*C1-cos(x)*log((1+sin(x))/cos(x))
 
 
 
Simplified =
 
sin(x)*C2+cos(x)*C1-cos(x)*log((1+sin(x))/cos(x))
 
 
 
                                                   1 + sin(x)
                sin(x) C2 + cos(x) C1 - cos(x) log(----------)
                                                     cos(x)

ODE =

x*D2y + Dy = x+1

 
Solution =
 
1/4*x^2+C1*log(x)+x+C2
 
 
 
Simplified =
 
1/4*x^2+C1*log(x)+x+C2
 
 
 
                               2
                          1/4 x  + C1 log(x) + x + C2

ODE =

x*D2y - 2*x*Dy + 2*y = x^2+2

 
Solution =
 
(-1/x*exp(2*x)-2*Ei(1,-2*x))*x*C2+x*C1+(Int((exp(2*x)+2*Ei(1,-2*x)*x)*(x^2+2)*exp(-2*x)/x,x)*x*exp(2*x)+(Ei(1,-2*x)*x+1/2*exp(2*x))*(x^2+x+5/2))*exp(-2*x)
 
 
 
Simplified =
 
-C2*exp(2*x)-2*C2*Ei(1,-2*x)*x+x*C1+1/2*x^2+exp(-2*x)*x^3*Ei(1,-2*x)+1/2*x+exp(-2*x)*Ei(1,-2*x)*x^2+5/4+5/2*exp(-2*x)*Ei(1,-2*x)*x+Int((exp(2*x)+2*Ei(1,-2*x)*x)*(x^2+2)*exp(-2*x)/x,x)*x
 
 
 
                                             /
  /  exp(2 x)                \               |
  |- -------- - 2 Ei(1, -2 x)| x C2 + x C1 + |
  \     x                    /               |
                                             \

          /                                2
         |  (exp(2 x) + 2 Ei(1, -2 x) x) (x  + 2) exp(-2 x)
         |  ----------------------------------------------- dx x exp(2 x)
         |                         x
        /

                                                        \
                                            2           |
         + (Ei(1, -2 x) x + 1/2 exp(2 x)) (x  + x + 5/2)| exp(-2 x)
                                                        |
                                                        /

ODE1 =

D2x - x - 2*y = t


ODE2 =

D2y - 2*y - 3*x = 1

 
Solution1 =
 
-C1*sin(t)-C2*cos(t)+2/3*C3*exp(2*t)+2/3*C4*exp(-2*t)-1/2+1/2*t
 
 
 
Solution2 =
 
1/4-3/4*t+C1*sin(t)+C2*cos(t)+C3*exp(2*t)+C4*exp(-2*t)
 
 
 
Simplified1 =
 
-C1*sin(t)-C2*cos(t)+2/3*C3*exp(2*t)+2/3*C4*exp(-2*t)-1/2+1/2*t
 
 
 
Simplified2 =
 
1/4-3/4*t+C1*sin(t)+C2*cos(t)+C3*exp(2*t)+C4*exp(-2*t)
 
 
 
   -C1 sin(t) - C2 cos(t) + 2/3 C3 exp(2 t) + 2/3 C4 exp(-2 t) - 1/2 + 1/2 t
 
       1/4 - 3/4 t + C1 sin(t) + C2 cos(t) + C3 exp(2 t) + C4 exp(-2 t)

ODE1 =

Dx - y = t


ODE2 =

Dy - x = 1

 
Solution1 =
 
exp(t)*C2-exp(-t)*C1-2
 
 
 
Solution2 =
 
exp(t)*C2+exp(-t)*C1-t
 
 
 
Simplified1 =
 
exp(t)*C2-exp(-t)*C1-2
 
 
 
Simplified2 =
 
exp(t)*C2+exp(-t)*C1-t
 
 
 
                          exp(t) C2 - exp(-t) C1 - 2
 
                          exp(t) C2 + exp(-t) C1 - t

ODE1 =

Dx - y = t


ODE2 =

Dy - x = 1

 
Solution1 =
 
11/10*exp(t)+exp(-t)-2
 
 
 
Solution2 =
 
11/10*exp(t)-exp(-t)-t