Differential equations Many laws of science and engineering are expressed in differential equations (DE). A DE is an equation that involves one or more derivatives of an unknown function. Solving a DE means finding a solution that satisfies the DE. But only some simple differential equations can be solved by explicit formulas.

Classification differential equations Differential equations are classified in several ways. An ordinary differential equation (ODE) is one that involves derivatives with respect to only one variable. A partial differential equation (PDE) is one that involves partial derivatives of unknown functions with respect to more than one variable.

Differential equations are also classified according to the highest order of derivatives in an equation. If the highest order is 1, the DE is called first-order DE. A DE of second or higher order are called higher-order DE.

A DE of order \(n\) has a general solution with \(n\) constants. In general, the n-th order ODE can be written as \(F(x,y,y',y'',\cdots,y^{(n)})=0\), where \(x\) is the independent variable, and \(y\) is the unknown function. In particular, the first-order ODE can be expressed as \(F(x,y,y')=0\), and the second-order ODE as \(F(x,y,y',y'')=0\).

Examples A simple first-order ODE is of the form \(y'=g(x)\), where the antiderivative \(y=∫g(x)dx\) is the solution (if it exists) to the DE. For example, \(y'=3x^2, y'=2x^{-2}+x\) and \(y'=\cos2x-e^x+3\) are simple first-order DEs. Check itg;3*x**2;x || itg;2*x**(-2)+x;x || itg;cos(2*x)+exp(x)+3;x ||.

An equation of the form \(y'+p(x)y=q(x)\) is a standard form first-order linear differential equation, where \(p(x),q(x)\) are continuous functions of \(x\) on some common interval. They are also called non-constant coefficients.

A linear ODE is homogeneous if \(q(x)=0\); a linear ODE is non-homogeneous if \(q(x)≠0\). In the following sections, we focus on solving some special form of first- and second-order ODE.

Solutions to differential equations A function \(y=f(x)\) is a solution to a DE if the DE holds when \(y\) and its derivatives are replaced by \(f(x)\) and the derivatives of \(f(x)\).

Examples The exponential function \(y=e^{kt}\) satisfies the DE \(y'=ky\), so \(y=e^{kt}\) is called a particular solution to the DE. All functions of the form \(y=Ce^{kt}\) for any constant \(C\) not only satisfy the DE, but also contain all possible solutions to the DE, so \(y=Ce^{kt}\) is called the general solution to the DE. Check dif;C*exp(k*t);t || ode;y_1-k*y;x ||.

In many applications, we are interested of finding a particular solution to a DE that satisfies a set of initial conditions. The problem of finding a particular solution of a DE that satisfies some initial conditions is called initial-value problem. The number of initial conditions must match the number of constants in the general solution of a DE.

Ordinary differential equations: [ode; expr; iv], [ ode(expr, iv) ]

An ordinary DE contains an unknown function y(x) of one variable x, its derivatives y', y'' or other higher-order derivatives, and some functions of x such as p(x) and q(x). The unknown function y depends on x, and thus x is an independent variable "iv".

Use the code pattern "ode; expr; iv" to help find solutions to ODE, where "expr" is the expression of a differential equation excluding the portion "= 0", and "iv" is the independent variable involved in the ODE.

An ODE must include derivatives of an unknown function with respect to its independent variable, and you must write a valid ODE expression so that the program can take it. To write the derivatives of an unknown function y with respect to x, for example, write an underscore after y, followed by an integer indicating the order of derivatives. In this way, "y_1" stands for first derivative y', "y_2" for the second derivative y'', and "y_n" for the nth-order derivative y⁽ⁿ⁾ . Write the derivatives of unknown functions with other names in a similar fashion.

The module "ode" help find solutions to some forms of ODEs such as linear or nonlinear, homogeneous or heterogeneous, constant- or non-constant-coefficients, but second or higher order nonlinear ODEs are not supported.

Attention Use the code pattern "ode(expr, iv)" to get the same results as "ode; expr; iv".

Attention If a DE is expressed in terms of f(x), g(y), h(z) and their derivatives, you need to rewrite it using a single variable for functions and apply the above code patterns to find its solution. For instance, a DE g(x) + 3g'(x) = 2 can be written as y + 3y_1 = 2 or g + 3g_1 = 2, and then apply ode;y+3*y_1-2;x || ode(g+3*g_1-2,x) ||.

Attention Use 'dif" module to write the above DE g(x) + 3g'(x) = 2 as g(x) + 3*dif(g(x),x) - 2, and then apply ode(g(x)+3*dif(g(x),x)-2) || to get the same solution for g(x). Combining "ode" and "dif" modules for solving DE is conceptually straightforward. You must write the unknown function as g(x) and write variable x as a single letter. Don't mix using derivatives such as "g_1" and "dif(g(x),x)" in the same DE for g'(x). The DE written in this way clearly specifies each term of derivatives, unknown functions, and variables, making it ready for using the "ode" module.

Examples Solve for h(z) for h''(z) + 9h(z) = 0. Let y = h(z). Then we have y_2 + 9y = 0. Solve y by ode;y_2+9*y;z || ode(y_2+9*y,z) || ode;h_2+9*h;z || or ode(dif(h(z),z,2)+9*h(z)) ||.

Examples
(1) Solve \(y'-2x=0\) by ode;y_1-2*x;x || ode(y_1-2*x,x) || ode(dif(y(x),x)-2*x) ||.
(2) Solve \(g''(z)-g'(z)=z\) by ode;g_2-g_1-z;z || ode(g_2-g_1-z,z) || ode(dif(g(z),z,2)-dif(g(z),z)-z) ||.
(3) Solve \( y''-y'-2y=0\) by ode;y_2-y_1-2*y;x || ode(y_2-y_1-2*y,x) || ode(dif(y(x),x,2)-dif(y(x),x)-2*y(x)) ||.
(4) Solve \(y''+2y'+3y=\sin x\) by ode;y_2+2*y_1+3*y-sin(x);x || ode(y_2+2*y_1+3*y-sin(x),x) || ode(dif(y(x),x,2)+2*dif(y(x),x)+3*y(x)-sin(x)) ||.

Examples ode;y_1-x/y;x || ode;s_2-g;g || ode;y_1*y-1;x || ode;y_1+2*y-exp(-3*x);x || ode;2*y_1-x*y-2*x;x || ode;y_1-y/x-y**3;x || ode;z_1+3*z-exp(2*x);x || ode;v*u_1+u-cos(v);v || ode;y_1-x*exp(x);x || ode;y_1-y*exp(x);x || ode;y**3*(1+x)-y_1;x || ode;y_1-2*y-y**2;x || ode;y_2-y-x;x || ode;y_1+y-x*y**2;x || ode;y_4-2*y_2+y;x || ode;y_1-(x**2-3*x+2)/(2*y+5);x || ode;y_1+2*y-2*y**(1/2);x || ode;x*y_1-y-x*log(x);x || ode;x_2+2*x_1+6*x;t || ode;x_2+k*x;t || ode;y_2+6*y;x || ode;x_2-x_1+t**2;t || ode;y_2+4*a*x;x || ode;t**2*x_2-t*x_1-3*x;t || ode;y_2+y_1+y-exp(2*x);x || ode;y_2+2*y_1+2*y-x**3+cos(x);x || ode;y_2+2*y_1+3*y-exp(x)-x;x || ode;y_2-y_1-sin(x)+x;x || ode;y_2-2*y_1+5*y;x || ode;y_2-2*y_1+y-x;x || ode;y_2-5*y_1+3*y-cos(x);x || ode;y_2+4*y_1-y-exp(-2*x);x || ode;y_3+y_2-y_1-y;x || ode;x**2*y_4+1;x || ode;y_3-4*y_2+3*y_1;x || ode;x_4+2*x_2+x;t || ode;x_2-x_1/t+x/t**2;t || ode;x**2*y_2+x*y_1-y-x;x || ode;x**2*y_2+x*y_1-y-1/x;x || ode;x**2*y_2+x*y_1-y-x**r;x || ode;y_2-2*y_1+y;x || ode;y_2+3*y_1-2*y-x;x||.

Examples
(1) If an ODE is of the form \(ydx+xdy=0\), rewrite it as \(y+xy'=0\), and type ode;y+x*y_1;x || for the solution.
(2) For \(xdx+dy=e^xdx\), type ode;x+y_1-exp(x);x ||.
(3) For \(xdx+(1+x^2)dy=0\), type ode;x+(1+x**2)*y_1;x ||.
(4) For \((x+y)dx=(y-x)dy\), and type ode;x+y+(x-y)*y_1;x ||.
(5) For \(\frac{dx}{y}+\frac{2dy}{x}=0\), type ode;1/y+2*y_1/x;x ||.

Exponential functions
(1) The exponential function \(y=f(x)=e^x\) can be defined by the differential equation \(y'=y\) and an initial value \(f(0)=1\) for all \(x\). Check ode;y_1-y;x ||, which returns \(y=f(x)=Ce^x\), and \(f(0)=1⇒C=1\).

(2) The function \(y=f(x)=e^{-x}\) can be defined as \(y'=-y\) and \(f(0)=1\) for all \(x\). Check ode;y_1+y;x ||, which returns \(y=Ce^{-x}\), and \(f(0)=1⇒C=1\). In both cases, \(f(x) > 0\) and \(f'(x) > 0\) for all \(x\).

In this section, we focus first-order ODE on (1) separation of variables, (2) linear and homogeneous equations, (3) linear and non-homogeneous equations, (4) exact equations, and (5) Bernoulli equations.

Separation of variables If the function \(f(x,y)\) of the DE \(y'=f(x,y)=g(x)h(y)\) can be factored as a product of a function of \(x\) and a function of \(y\), the DE is called separable. Rewrite the DE as \(\frac{dy}{dx}=g(x)h(y)\) and \(\frac{dy}{h(y)}=g(x)dx\). Then integrate both sides \(\int \frac{dy}{h(y)}=\int g(x)dx\), complete the integration, and we can obtain the solution of \(y\).

In case the DE is of the form \(p(x)+q(y)y'=0⇒∫p(x)dx+∫q(y)dy=C⇒P(x)+Q(y)=C\), where \(P(x)\) is antiderivative of \(p(x)\) and \(Q(y)\) is antiderivative of \(q(y)\). Thus, the implicit equation is the solution to the DE.

Examples
(1) Solve \(\frac{dy}{dx}=2x+3xy⇒\frac{dy}{2+3y}=xdx⇒∫\frac{dy}{2+3y}=∫xdx⇒\frac{1}{3}\ln|2+3y|=\frac{x^2}{2}+C⇒\)\(y=\frac{C}{3}e^{\frac{3}{2}x^2}-\frac{2}{3}\). Check ode;y_1-x*(2+3*y);x ||.

(2) Solve \(y'y=2x⇒ydy=2xdx⇒∫ydy=∫2xdx\) \(⇒\frac{y^2}{2}=x^2+C⇒y^2-2x^2=2C\). Thus, the solution \(y\) is implicitly defined in the equation. Check ode;y_1*y-2*x;x ||.

(3) Solve \(y'=\frac{2y-1}{3x+4}=0⇒\frac{dy}{2y-1}=\frac{dx}{3x+4}⇒∫\frac{dy}{2y-1}=∫\frac{dx}{3x+4}⇒\frac{1}{2}\ln|2y-1|= \frac{1}{3}\ln|3x+4|+C⇒\) \(\ln\frac{|2y-1|} {\sqrt[3]{(3x+4)^2}}=2C⇒|2y-1|=e^{2C}(3x+4)^{\frac{2}{3}}⇒y=C_1(3x+4)^{\frac{2}{3}}+\frac{1}{2}\) by ode;y_1-(2*y-1)/(3*x+4);x ||.

Examples ode;y_1-3*y/(2*x);x || ode;x_1-t**2/x**2;t || ode;y_1-x**2*y**2;x || ode;y_1-exp(y)*cos(x);x || ode;x_1-1-x**2;t || ode;y_1-t*(1+exp(y));x || ode;y_1-6*x**2*y;x || ode;y_1-exp(-y)*sin(x);x ||.

First-Order Homogeneous ODE

Homogeneous functions A function \(f(x,y)\) is homogeneous of degree \(n\) if \(f(tx,ty)=t^nf(x,y)\) for some constant \(n\) and real number \(t\). For example, \(f(x,y)=2x^2+3xy+4y^2\) is homogeneous of degree 2 because \(f(tx,ty)=t^2f(x,y)\). The function \(f(x,y)=\frac{y}{x}\) is homogeneous of degree 0. But \(f(x,y)=3x + y^3\) is not homogeneous.

Homogeneous differential equations A DE \(M(x,y)dx+N(x,y)dy=0\) is called homogeneous if \(M(x,y)\) and \(N(x,y)\) are homogeneous functions of the same degree.

Examples
(1) The DE \((2x^2y+3x^3)dx+2(y^2x+y^3)dy=0\) is homogeneous.
(2) The DE of the form \(y'=f(\frac{y}{x})\) is a homogeneous DE.
(3) \(e^{\frac{y}{x}}dx+\sin\frac{y}{x}dy=0\) is not a homogeneous DE.

We can solve a homogeneous DE of the form \(y'=f(\frac{y}{x})\) by substituting \(v(x)=\frac{y}{x}\), where \(v(x)\) is a differentiable function of \(x\). Differentiating both sides of \(y=v(x)x\), we get \(y'=v'(x)x+v(x)\). Substituting \(y,y'\) in the DE \(y'=f(\frac{y}{x})⇒v+v'x=f(v)⇒v'x=f(v)-v⇒\frac{dv}{f(v)-v}=\frac{dx}{x}\) by separation of variables.

Examples
(1) Solve \(y'=\frac{y^2}{xy+x^2}\). Since \(\frac{y^2}{xy+x^2}=\frac{(y/x)^2}{y/x+1}=\frac{v^2}{v+1},\frac{dv}{f(v)-v}= -\frac{v+1}{v}dv=\frac{dx}{x}⇒-∫\frac{v+1}{v}dv=∫\frac{dx}{x}⇒\)
\(-(v+\ln|v|)=\ln|x|+C⇒\ln|x|+\ln|v|=-C-v⇒\ln|y|=-C-v⇒y=±e^{-C-\frac{y}{x}}⇒y=C_1e^{-\frac{y}{x}}\). Check ode;y_1-y**2/(x*y+x**2);x ||.

(2) Solve \((x^2-y^2)dx+2xydy=0\). It is homogeneous but hard to solve by separation of variables. Let \(y=v(x)x\). Then \(dy=xdv+vdx\). Substituting \(y,dy\) into the DE, we have \((x^2v-v^2x^2)dx+\)\(2vx^2(xdv+vdx)=0\). Simplifying the equation, we get \(\frac{dx}{x}=\frac{-2v}{1+v^2}dv\). Integrating both sides, we obtain \(\ln|x|=-\ln(1+v^2)+C\), or \(\ln|x|(1+v^2)=C⇒|x|(1+\frac{y^2}{x^2})= e^C⇒x^2+y^2=C_1x\) by ode;x**2-y**2+2*x*y*y_1;x ||.

First-order linear homogeneous ODE An equation of the form \(y'+p(x)y=q(x)\) is called a first order linear DE in standard form, where \(p(x), q(x)\) are continuous functions of \(x\) on some interval.

If \(p(x)=0\), the equation reduces to \(y'=q(x)\), and the solution of this DE is the antiderivative \(y=∫q(x)dx\).

If \(q(x)=0\), the DE \(y'+p(x)y=0\) is called linear homogeneous DE. By separation of variables, \(\frac{y'}{y}=-p(x)⇒\)\(∫\frac{dy}{y}=∫-p(x)dx⇒\ln|y|=-∫p(x)dx\). So the general solution of the equation is \(y = Ce^{-\int p(x)dx}\).

In particular if \(p(x)=k\) is a constant, the DE becomes \(y'+ky=0\), and the general solution to the DE is \(y = Ce^{-kx}\). Solving the initial value problem \(y'=-ky, y(0)=C\), we have the unique solution \(y=y_0e^{-kt}\). The DE \(y'=ky\) is often used for exponential growth modeling.

Examples
(1) Solve \(y'-2xy=0\). Since \(p(x)=-2x,v(x)=∫2xdx=x^2\) is the integrating factor, and the general solution is \(y=Ce^{x^2}\). Check ode;y_1-2*x*y;x ||.

(2) Solve \(xy'-y=0\). Convert the DE to standard form \(y'-\frac{y}{x}=0,p(x)=-\frac{1}{x}\), and the integrating factor is \(v(x)=∫\frac{1}{x}dx=\ln|x|\). So the general solution is \(y=Ce^{\ln|x|}=C|x|=C_1x\) for constant \(C\). || ode;x*y_1-y;x ||.

(3) Solve \(y'\cos x+y\sin x=0,⇒y'+y\tan x=0\)\(⇒y=Ce^{-∫\tan x dx}=Ce^{\ln |\cos x|}=C_1\cos x\). Check ode;y_1+tan(x)*y;x || ode;cos(x)*y_1+sin(x)*y;x ||.

First-Order Non-homogeneous ODE

First-order linear non-homogeneous ODE If \(q(x)≠0\), the standard linear DE \(y'+p(x)y=q(x)\) is non-homogeneous. We use two methods to solve the non-homogeneous DE.

(1) Variation of parameters Given \(y = Ce^{-\int p(x)dx}\) is the solution to the homogeneous DE \(y'+p(x)y=0\), we assume \(y=C(x)e^{-\int p(x)dx}\) is the solution to the non-homogeneous DE \(y'+p(x)y=q(x)\). Now we need to determine the function \(C(x)\).

Since \(y'=e^{-\int p(x)dx}[C'(x)-C(x)p(x)]\), substituting \(y,y'\) in the non-homogeneous DE, we have \(e^{-\int p(x)dx}[C'(x)-C(x)p(x)]+p(x)C(x)e^{-\int p(x)dx}=q(x)⇒C'(x)=q(x)e^{\int p(x)dx}⇒\)
\(C(x)=∫q(x)e^{\int p(x)dx}dx⇒y=e^{-\int p(x)dx}∫q(x)e^{\int p(x)dx}dx\) is the general solution to the non-homogeneous DE.

(2) Integrating factor We can think of the left side of the DE is the derivative \([v(x)·y]'\), the product of \(y\) and some function of \(x\). Suppose there exists such an integrating factor \(v(x)\ne 0\) and \(v(x)\) is a differentiable function of \(x\) in the domain of \(y(x)\). Then \([v(x)y]'=v(x)q(x)\). Solving \([v(x)y]'=v(x)q(x)\) by integrating both sides, we obtain \(v(x)y=\int v(x)q(x)dx\) or \(y=\frac{1}{v(x)}\int v(x)q(x)dx\), which is the solution of the DE in terms of \(v(x)\) and \(q(x)\).

To find \(v(x)\), calculate and match the derivatives \([v(x)y]'=v'(x)y+v(x)y'=v(x)[y'+p(x)y]=v(x)y'+v(x)p(x)y\). Simplifying the equation, we have \(v'(x)=v(x)p(x)\) or \(\frac{v'(x)}{v(x)}=p(x)\). Integrating both sides, we obtain \(\int \frac{v'(x)}{v(x)}dx=\int p(x)dx\) and \(\ln|v(x)|=\int p(x)dx\), or \(v(x)=\pm e^Ce^{\int p(x)dx}\).

We don't need a general expression for the factor. So we take a simple positive integrating factor for \(v(x)=e^{\int p(x)dx}\). Substituting \(v(x)\) into the expression of \(y\), we find the solution of the differential equation by \(y=e^{-\int p(x)dx}\int q(x)e^{\int p(x)dx}dx\).

Do not memorize these formulas. The key is to find the integrating factor. Memorize its form and use it to complete the integration for the final solution. If an equation is not in standard form, we need to convert it to a standard form and follow the above procedure to find the solution.

Examples
(1) Solve \(y'-2xy=3x\). Since \(p(x)=-2x,q(x)=3x\), the general solution is \(y=e^{∫2xdx}∫3x e^{∫-2xdx}\) \(=e^{x^2}∫3xe^{-x^2}dx=Ce^{x^2}-\frac{3}{2}\). Check ode;y_1-2*x*y-3*x;x ||.

(2) Solve \(xy'+x^2y=2x^2\). Convert the equation to standard form \(y'+xy=2x\), whose solution is \(y=e^{-∫xdx}∫2xe^{∫xdx}dx=e^{-\frac{x^2}{2}}∫2xe^{\frac{x^2}{2}}dx=e^{-\frac{x^2}{2}}(2e^{\frac{x^2}{2}}+C)\)\(=2+Ce^{-\frac{x^2}{2}}\) by ode;x*y_1+x^2*y-2*x^2;x ||.

(3) Solve \(xy'+y=e^x\). Convert it to a standard form \(y'+\frac{y}{x}=\frac{e^x}{x}\). The general solution is \(y=e^{∫\frac{-1}{x}dx}∫\frac{e^x}{x}e^{∫\frac{1}{x}dx}dx=e^{-\ln|x|}∫\frac{e^x}{x}e^{\ln|x|}dx=\frac{1}{x}(e^x+C)\) for \(x≠0\). Check ode;x*y_1+y-exp(x);x ||.

First-order exact ODE A DE of the form \(M(x,y)dx+N(x,y)dy=0\) or \(\frac{dy}{dx}=-\frac{M(x,y)}{N(x,y)}\) is called an exact DE if there is a function \(z=f(x,y)\) whose total differential \(dz=\frac{∂f}{∂x}dx+\frac{∂f}{∂y}dy=M(x,y)dx+N(x,y)dy\). The level curves \(f(x,y)=C\) are the solution curves to the DE.

Recall the total differential of a two-variable function \(z=f(x,y)\) is \(dz=\frac{∂f}{∂x}dx+\frac{∂f}{∂y}dy=f_xdx+f_ydy\). In other words, \(M(x,y)dx+N(x,y)dy=0\) is an exact DE if there is a potential function \(f\) whose gradient is \(Δf=⟨M(x,y),N(x,y)⟩\), or if the vector field \(M(x,y){\bf i}+N(x,y){\bf j}\) is conservative. Now solving the exact DE is equivalent to solving \(dz=0\) or \(Δf=0\) for a potential function \(f\) of the conservative field \(⟨M(x,y),N(x,y)⟩\). Since \(Δf=0⇒f(x,y)=C\) for some constant \(C\).

Note that the field \(⟨M(x,y),N(x,y)⟩\) is conservative if and only if their cross-partial derivatives are equivalent \(\frac{∂M}{∂y}=\frac{∂N}{∂x}\). Thus, a DE of the form \(M(x,y)dx+N(x,y)dy=0\) is an exact DE at every point \((x,y)\) in a simply connect open domain if and only if \(\frac{∂M}{∂y}=\frac{∂N}{∂x}\). We use this fact to determine if a DE is an exact DE.

Example Solve \((x+3x^2y)dx+(x^3+3y^2)dy=0\). Since \(M(x,y)=x+3x^2y,N(x,y)=x^3+3y^2\), we have \(\frac{∂M}{∂y}=3x^2,\frac{∂N}{∂x}=3x^2\). So the DE is an exact DE, and there is \(z=f(x,y)\) such that \(f_x=M,f_y=N\). Check the vector csv;(x+3*x**2*y)*i+(x**3+3*y**2)*j;x;y || cul;(x+3*x**2*y)*i+(x**3+3*y**2)*j;x;y ||.

Since \(f_x=x+3x^2y⇒f(x,y)=∫(x+3x^2y)dx=\frac{x^2}{2}+x^3y+g(y)\). To determine \(g(y)\), we have \(f_y=x^3+g'(y)=N=x^3+3y^2⇒\)\(g'(y)=3y^2⇒g(y) =∫3y^2dy=y^3+C\). Thus, the potential function or solution is \(f(x,y)=\frac{x^2}{2}+x^3y+y^3+C=0\). The general solution to the DE is the implicit equation \(\frac{x^2}{2}+x^3y+y^3+C=0\) for some constant \(C\). Check ode;x+3*x**2+(x**3+3*y**2)*y_1;x ||, which gives a form solution that is hard to compare with the implicit equation. However, we can verify that \(\frac{dy}{dx}=-\frac{M(x,y)}{N(x,y}=\frac{-(x+3x^2y)}{x^3+3y^2}\) by the implicit differentiation dif;x**2/2+x**3*y+y**3+C;x;y ||.

First-order Bernoulli ODE A first-order nonlinear DE of the form \(y'+p(x)y=q(x)y^n\) is called Bernoulli equation for a real number \(n≠0, n≠1\). The DE is linear for \(n=0\), and the DE is separable for \(n=1\).

To find the solution, multiplying both sides of the Bernoulli equation by a term \((1-n)y^{-n}\), we have the following equation \((1-n)y^{-n}y'+(1-n)p(x)y^{1-n}=(1-n)q(x)\). Notice that the first term on the left side is equal to \(\frac{d}{dx}y^{1-n}\) and the DE becomes a first-order linear DE.

To simplify the notation, let \(z=y^{1-n}\). Then \(z'+(1-n)p(x)z=(1-n)q(x)\). So the integrating factor for this linear DE is \(v(x)=e^{-∫(1-n)p(x)}\), and the solution to the linear DE is \(z=\frac{1}{v(x)}\int v(x)(1-n)q(x)dx\).

Substituting \(v(x)\) and \(z\) in the solution, we have \(y^{1-n}=e^{-(1-n)\int p(x)dx}\int (1-n)q(x)e^{(1-n)\int p(x)dx}dx\). Completing the integration, we can obtain the solution of \(y\) to the Bernoulli equation.

Example
(1) Solve \(y'+2y=3e^xy^2\). Since \(p(x)=2,q(x)=3e^x,n=2\), the solution to the Bernoulli DE is \(y^{-1}=e^{∫2dx}∫-3e^xe^{∫-2dx}dx=e^{2x}(3e^{-x}+C)=3e^x+Ce^{2x}⇒y=\frac{1}{3e^x+Ce^{2x}}\), ode;y_1+2*y-3*exp(x)*y**2;x ||.

(2) Solve \(y'+xy=2xy^3\). Given \(p(x)=x,q(x)=2x,n=3\), the general solution satisfies \(y^{-2}=e^{∫2xdx}∫-4xe^{-∫2xdx}dx=e^{x^2} (2e^{-x^2}+C)=2+Ce^{x^2}⇒y^2=\frac{1}{2+Ce^{x^2}}\) by ode;y_1+x*y-2*x*y^3;x.

Example More Bernoulli equations by ode;x*y_1-y-x*y**2;x || ode;y_1-y-x*y**(-1);x || ode;x*y_1+y-y**2*log(x);x ||.

Second-order ODE We focus the second-order ODE on the following three forms of DE
(1) Some simple forms such as \(y''=f(x)\), \(y''=f(x,y')\) and \(y''=f(y,y')\).
(2) Linear constant-coefficient homogeneous equations of the form \(y''+by'+cy=0\).
(3) Linear constant-coefficients non-homogeneous equations \(y''+by'+cy=f(x)\), where \(f(x)\) has some special forms of polynomials, since and cosine, and exponential functions, or a simple combination of these functions.

We first introduce four basic theorems regarding the solutions to the second order ODE, and then use these theorems to solve second-order linear constant coefficients homogeneous ODE and non-homogeneous ODE of some special forms.

Theorem 1 If \(y_1,y_2\) are solutions to the DE \(y''+p(x)y'+q(x)y=0\), then any linear combination of \(y_1\) and \(y_2\) is also a solution to the DE.

To see this, let \(y_p=Ay_1+By_2\) for any constants \(A,B\). Then \(y_p''=Ay_1''+By_2'',y_p'=Ay_1'+By_2'\) and \(y_p''+p(x)y_p'+q(x)y_p=A(y_1''+p(x)y_1'+q(x)y_1)+B(y_2''+p(x)y_2'+q(x)y_2)=A(0)+B(0)\) = 0.

Theorem 2 If \(y_1,y_2\) are (linearly independent) solutions to the DE \(y''+p(x)y'+q(x)y=0\), and \(y_1(x)≠ky_2(x)\) for any constant \(k\), then the linear combination \(Ay_1+By_2\) is the general solution to the DE.

Example \(y_1=\sin x,y_2=\cos x\) are solutions to \(y''+y=0\), so \(A\sin x+B\cos x\) is the general solution to the DE \(y''+y=0\). In particular, \(\sin(x-2)=\sin x\cos2-\cos x\sin2\) and \(\cos(x+3)=\cos x\cos3-\sin x\sin3\) both are solutions to the DE, but \(\sin2x\) and \(\cos 2x\) are not solutions to the DE.

Theorem 3 If \(y_p\) is a particular solution to the DE \(y''+p(x)y'+q(x)y=f(x)\) and \(Y\) is the general solution to the DE \(y''+p(x)y'+q(x)y=0\), then the sum \(y_p+Y\) is the general solution to the non-homogeneous DE.

Substituting \(y_p+Y,\quad y'_p+Y'\), and \(y''_p+Y''\) in the non-homogeneous DE, we obtain \(y_p''+Y''+p(x)(y_p'+Y')\)
\(+q(x)(y_p+Y)=y_p''+p(x)y_p'+q(x)y_p+Y''+p(x)Y'+q(x)Y=f(x)+0=f(x)\). Since the solution \(y=y_p+Y=y_p+Ay_1+By_2\) has two constants, \(y\) is the general solution.

Example \(y_p=x^2\) is a particular solution to \(y''+y=x^2+2\), the general solution is \(y=A\cos x+B\sin x+x^2\).

Theorem 4 If \(y_1,y_2\) are particular solutions to \(y''+p(x)y'+q(x)y=f(x)\) and \(y''+p(x)y'+q(x)y=g(x)\), respectively, then \(y_1+y_2\) is a particular solution to \(y''+p(x)y'+q(x)y=f(x)+g(x)\). This theorem is called superposition principle.

Examples
(1) Solve \(y''=f(x)\). Let \(y'=p(x)\). Then \(y'=∫f(x)dx+C,y=∫p(x)dx=∫[∫f(x)dx+C]dx\). Check ode;y_2+cos(x)-1;x || for \(y''=1-\cos x⇒y'=x-\sin x+C⇒\)\(y=\frac{x^2}{2}+\cos x+Cx+C_1\).

(2) Use a similar method to solve higher-order DE of the form \(y^{(n)}=f(x)\). Let \(y'''=2x\). Then \(y'''=2x⇒y''=x^2+C⇒y'=\frac{x^3}{3}+Cx+C_1\)\(⇒y=\frac{x^4}{12}+\frac{Cx^2}{2}+C_1x+C_2\). Check ode;y_3-2*x;x ||.

(3) Solve \(y''=f(x,y')\). Let \(y'=p(x)\). Then \(y''=p'=f(x,p)⇒p=y'=∫f(x,p)dx=g(x,C)\)
\(⇒y= ∫g(x,C)dx+C_1\). Check ode;y_2-y_1+2*x;x || for \(y''=y'-2x\). Let \(y'=p,y''=p'\). Then \(p'-p=-2x⇒p=e^{∫1dx} ∫-2xe^{∫-1dx}dx=e^x(∫-2xe^{-x}dx)=e^x[2e^{-x}(x+1)+C]\)
\(=2(x+1)+Ce^x⇒y=∫[(2x+1)+Ce^x]dx=x^2+2x+Ce^x+C_1\).

(4) Solve \(y''=f(y,y')\). Let \(y'=p\). Then \(y''=p'=\frac{dp}{dy}\frac{dy}{dx}=\frac{dp}{dy}p=f(y,p)⇒p=\frac{dy}{dx}=∫f(y,p)dy\)
\(=g(y,C)⇒\frac{dy}{g(y,C)}=dx⇒x=∫\frac{dy}{g(y,C)}+C_1\).

Second-order linear constant homogeneous ODE The DE of the form \(ay''+by'+cy=0\) is a second-order, linear homogeneous DE with some constant coefficients \(a≠0,b,c\).

Just as the first-order constant coefficient homogeneous DE \(y'-ky=0\) having solution \(y=Ce^{kt}\), we also assume \(y=e^{kt}\) is a solution to the DE \(ay''+by'+cy=0\). Substituting it to the DE, we have \(e^{kt}(ak^2+bk+c)=0\). Since \(e^{kt}\) will never be 0, the equation \(ak^2+bk+c=0\), which called characteristic (auxiliary) equation, must have solutions.

If the \(e^{kt}\) is a solution to the DE, \(k=\frac{-b\pm\sqrt{b^2-4ac}}{2a}=\frac{-b\pm\sqrt{D}}{2a}\) must be the root(s) to the characteristic equation. So solving the DE \(ay''+by'+cy=0\) boils down to finding roots of the quadratic equation \(ak^2+bk+c=0\).

Case 1 If \(D>0\), then \(k_1=\frac{-b+\sqrt{D}}{2a}, k_2=\frac{-b-\sqrt{D}}{2a}\). Then the general solution is \(y=Ae^{k_1t}+Be^{k_2t}\).

Case 2 If \(D=0, k_1=k_2=-\frac{b}{2a}=k\). The solution is \(y=Ae^{kt}+Bte^{kt}\). Since \(y_1(t)=e^{kt}\) is one solution, the other can be written as \(y_2(t)=y_1(t)u(t)=e^{kt}u(t)\). Substituting \(y_2',y_2''\) into the DE and solving for \(u''(t)=0⇒u(t)=t\), we have \(y_2(t)=te^{kt}\), and the general solution is \(y=Ae^{kt}+Bte^{kt}\).

Case 3 If \(D< 0, k=-\frac{b}{2a}, β=\frac{\sqrt{D}}{2a}\), and the roots are \(r=k\pm i β\). So the two particular solutions to the DE are \(y_1(t)=e^{(k+i β)t},y_2(t)=e^{(k-i β)t}\), but they are not real valued. Using Eula's formula, \(y_1=e^{kt}(\cosβt+i\sinβt)\), \(y_2=e^{kt}(\cos βt-i\sin βt)\). Let \(y_3=\frac{y_1+y_2}{2}=e^{kt}\cos βt,y_4=\frac{y_1-y_2}{2i}=e^{kt}\sin βt\). Since \(y_3\) and \(y_4\) are another two solutions and \(\frac{y_3}{y_4}\) is not a constant, we have the general solution \(y=Ae^{kt}\cos( β t)+Be^{kt}\sin( βt)\).

Example
(1) Solve \(y''+2y'-3y=0\). The characteristic equation is \(k^2+2k-3⇒k_1=-3,k_2=1\), and the general solution is \(y=Ae^{-3x}+Be^x\). If \(y(0)=1,y'(0)=5\) then \(A+B=1\) and \(-3A+B=5⇒\)\(A=-1,B=2\). So the particular solution is \(y=-e^{-3x}+2e^x\). Check ode;y_2+2*y_1-3*y;x ||.

(2) Solve \(y''-4y'+4y\). The characteristic equation is \(k^2-4k+4=0⇒k=2\), and the general solution is \(y=Ae^{2x}+Bxe^{2}x\). Check ode;y_2-4*y_1+4*y;x ||.

(3) Solve \(y''+2y'+2y=0\). The characteristic equation is \(k^2+2k+2=0⇒k=-1±i\), and the general solution is \(y=Ae^{-x}\cos x+Be^{-x}\sin x\). Check ode;y_2+2*y_1+2*y;x ||.

(4) Based on Newton's second law of motion and Hooke's law, the equation of simple harmonic motion of a mass \(m\) on a horizontal spring is \(F=m\frac{d^2y}{dt^2}=-ky\), where \(y\) is its displacement from the equilibrium position, and \(k\) is a spring constant. Solve \(y''+\frac{k}{m}y=0\). The characteristic equation \(r^2+wr=0⇒r=±i\sqrt{w}\), and the general solution to the DE is \(y=A\cos wt+B\sin wt\) for \(w=\sqrt{\frac{k}{m}}\). Check ode;y_2+k*y/m;x ||.

(5) ode;2*y_2-3*y_1+4*y;x || ode;y_2-3*y_1+7*y;x || ode;-y_2-3*y_1+y;x || ode;y_2-6*y_1+9*y;x || ode;y_2+3*y;x ||. ode;3*x_2-x;t ||.

Second-order linear constant non-homogeneous ODE A second order linear non-homogeneous DE is of the form \(y''+py'+q =f(x)\), which is in contrast to the corresponding homogeneous DE \(y''+py'+q =0\), for some constants \(p, q\). If \(y_1\) and \(y_2\) are solutions to the non-homogeneous DE, then \(y_1-y_2\) is the solution to the homogeneous DE, for \((y_1''+py_1'+qy_1)-(y_2''+py_2'+qy_2)=0\)
\(⇒(y_1''-y_2'')+p(y_1'-y_2')+q(y_1-y_2)\) = 0. We use this fact to find the general solution to the non-homogeneous DE.

Suppose \(y_p\) is one solution to the non-homogeneous DE and \(Y\) is the general solution to the corresponding homogeneous DE. Then \(y=Y+y_p\) is the general solution to the non-homogeneous DE \(y''+py'+q =f(x)\).

To see this, let \(y\) be another solution to the non-homogeneous DE, then \(y-y_p=Y\) is the solution to the homogeneous DE. Since \(Y\) has two constants, \(y\) is the general solution. Now the question becomes how to determine a particular solution \(y_p\) to the non-homogeneous DE.

The general solution \(y\) depends on the form of \(f(x)\) and the roots of the characteristic equation corresponding to the homogeneous DE. We discuss three different forms of \(f(x)\) by use of the method of undetermined coefficients: simple polynomials, sine or cosine, and exponential functions, or their simple combinations.

I. Simple polynomials If \(f(x)=P_n(x)\) is a polynomial and if both roots \(k_1k_2≠0\) of the characteristic equation corresponding to \(y''+py'+q =0\), then \(y_p=W_n(x)\).

If \(k_1=0\) or \(k_2=0\), then \(y_p=xW_n(x)\); if \(k_1 =k_2=0, y_p=x^2W_n(x)\).

Example
(1) For a simple case, it is known that \(y=x\) is a solution to \(y''+y=x\), and \(y=A\cos x+B\sin x\) is the general solution to \(y''+y=0\). Thus, the solution to \(y''+y=x\) is \(y=x+A\cos x+B\sin x\) by ode;y_2+y-x;x ||.

(2) The general solution to \(y''+y=c\) for some constant \(c\) is \(y=c+A\cos x+B\sin x\). Check ode;y_2+y-c;x ||.

(3) Following the same logic, the general solution \(y=Ae^{-x}+Be^{2x}-1\) to \(y''-y'-2y=2\) since \(y=-1\) is a particular solution to \(y''-y'-2y=2\) and \(y=Ae^{-x}+Be^{2x}\) is the general solution to \(y''-y'-2y=0\). Check ode;y_2-y_1-2*y-2;x ||.

Examples
(1) Solve \(y''+2y'-3y=x^2\). The characteristic equation \(k^2+2k-3=0⇒k_1=-3,k_2=1\). Let \(y_p=Ax^2+Bx+C\) with unknown coefficients. Using the method of undetermined coefficients to find \(A,B,C\). Substituting \(y_p'=2Ax+B,y_p''=2A\) in the equation, we have \(-3Ax^2+(4A-3B)x+(2A+2B-3C)=x^2⇒-3A=1,A+2B-3C=0,4A-3B=0⇒A=\frac{-1}{3},B=\frac{-4}{9},C=\frac{-14}{27}\). Thus, \(y_p=\frac{-x^2}{3}-\frac{4x}{9}- \frac{14}{27}\), and the solution is \(y=C_1e^{-3x}+C_2e^x-\frac{x^2}{3}-\frac{4x}{9}-\frac{14}{27}\) by ode;y_2+2*y_1-3*y-x**2;x ||.

(2) Solve \(y''-y'=2x-5\). Since \(k^2-k=0⇒k=1;k=0\). Let \(y_p=x(Ax+B)\). Then \(y_p'=2Ax+B,y_p''=2A\). Substituting them in the equation \(2A-2Ax-B=2x-5⇒-2A=2⇒A=-1\) and \(2A-B=-5⇒B=3\). So \(y_p=-x^2+3x\), and the general solution is \(y=C_1e^x+C_2-x^2+3x\). Check ode;y_2-y_1-2*x+5;x ||.

(3) Solve \(y''=3x^2+2x\). The roots of the characteristic equation is \(k^2=0⇒k_1=k_2=0\). Let \(y_p=x^2(Ax^2+Bx+C)⇒\)\(y_p'=4Ax^3+3Bx^2+2Cx⇒y_p''= 12Ax^2+6Bx+2C ⇒12A=3\),
\(6B=2,2C=0⇒A=\frac{1}{4},B=\frac{1}{3},C=0\). Thus, \(y_p=\frac{x^4}{4}+\frac{x^3}{3}\), and the general solution is \(y=C_1x+C_2+\frac{x^4}{4}+ \frac{x^3}{3}\). Check ode;y_2-3*x**2+2*x;x ||.

In fact, \(y'=∫(3x^2+2x)dx=x^3+x^2+C_1⇒y=∫(x^3+x^2+C_1)dx=\frac{x^4}{4}+\frac{x^3}{3}+C_1x+C_2\).

II. Exponential functions If \(f(x)=ae^{ax}\) and if \(k_1,k_2≠ a, y_p=Ae^{ax}\); if \(k_1=a\) or \(k_2=a\), then \(y_p=Axe^{ax}\); if \(k_1 =k_2=a\), then let \(y_p=Ax^2e^{ax}\).

Example
(1) Solve \(y''-3y'+2y=3e^{2x}\). The characteristic equation is \(k^2-3k+2=0⇒k=2\), \(k=1\). Let \(y_p=Axe^{2x}\). Then \(y_p'=Ae^{2x}+2Axe^{2x}\), \(y_p''=4Ae^{2x}+4Axe^{2x}\). Substituting them in the equation and matching coefficients, we have \(A=3,y_p=3xe^{2x}\), and the general solution is \(y=C_1e^{2x}+C_2e^x+3xe^{2x}\). Check the result by ode;y_2-3*y_1+2*y-3*exp(2*x);x ||.

(2) Solve \(y''-4y=2e^{x}\). The characteristic equation is \(k^2-4=0⇒k=±2\). Let \(y_p=Ae^x\). Then \(y_p'=Ae^x=y_p''\). Substituting them in the equation and matching coefficients, we have \(A=-\frac{2}{3}\), \(y_p=-\frac{2}{3}e^x\), and the general solution is \(y=C_1e^{2x}+C_2e^{-2x}-\frac{2}{3}e^x\). Check ode;y_2-4*y-2*exp(x);x ||.

(3) Solve \(y''-4y'+4y=5e^{2x}\). The characteristic equation is \(k^2-4k+4=0⇒k_1=k_2=2\). Let \(y_p=Ax^2e^{2x}\). Then \(y_p'=2Axe^{2x}+2Ax^2e^{2x},y_p''= 6Axe^{2x}+2Ae^{2x}\). Substituting them in the equation and matching coefficients, we have \(A=\frac{5}{2},y_p=\frac{5}{2}x^2e^{2x}\), and the general solution is \(y=C_1e^{2x}+C_2xe^{2x}+\frac{5}{2}x^2e^{2x}\). Check ode;y_2-4*y_1+4*y-5*exp(2*x);x ||.

III. Sine and cosine functions Let \(f(x)=a\sinβ x\) or \(f(x)=b\cosβ x\). Then we deal with two cases.
(1) If \(k_1,k_2≠0+βi,y_p=A\cos β x+B\sin β x\).
(2) If \(k_1,k_2= 0+βi\), \(y_p=x(A\cosβ x+B\sinβ x)\).

Example
(1) Solve \(y''+y=\cos x\). The characteristic equation is \(k^2+1=0⇒k=±i\). Let \(y_p=x(A\cos x\)\(+B\sin x)\). \(y_p'=A\cos x+B\sin x+x(-A\sin x+B\cos x),y_p''=2(-A\sin x+B\cos x)-x(A\cos x+B\sin x)\). Substituting \(y_p,y''_p\) in the DE and matching the coefficients, we have \(A=0,B=\frac{1}{2},y_p=\frac{x}{2}\sin x\), and the general solution is \(y=C_1\cos x+C_2\sin x+\frac{1}{2}x\sin x\). Check ode;y_2+y-cos(x);x ||.

(2) Solve \(y''+4y=3\sin x\). The characteristic equation is \(k^2+4=0⇒\) \(k=±2i\). Let a particular solution be \(y_p=A\cos x+B\sin x\). Then \(y_p'=-A\sin x+B\cos x,y_p''=-A\cos x-B\sin x\). Substituting \(y_p,y''_p\) in the DE and matching coefficients, we have \(A=0,B=1,y_p=\sin x\), and \(y=C_1\cos2x+C_2\sin2x+\sin x\) is the general solution to the DE . Check ode;y_2+4*y-2*sin(x);x ||.

(3) Solve \(y''-2y'+5y=13\cos 3x\). The characteristic equation is \(k^2-2k+5=0⇒k=1±2i\). Let one solution be \(y_p=A\cos3x+B\sin3x,y_p'=-3A\sin3x+3B\cos3x\), \(y_p''= -9A\cos3x-9B\sin 3x\). Substituting \(y_p,y'_p,y''_p\) in the DE and matching coefficients, we have \(A=-1,B=\frac{-3}{2}\), \(y_p=-\cos3x-\frac{3}{2}\sin3x\), and the general solution is \(y=C_1e^x\cos2x+C_2e^x\sin2x-\cos3x-\frac{3}{2}\sin3x\). Check ode;y_2-2*y_1+5*y-13*cos(3*x);x ||.

IV. Mixed types of functions In case \(y''+py'+qy=f(x)\) and \(f(x)=e^{ax}[P_m(x)\cos β x+Q_n(x)\sin β x]\), where \(P_m(x)\) and \(Q_n(x)\) are polynomials.
(1) If both roots \(k_1,k_2≠ a+β i\) of the characteristic equation corresponding to \(y''+py'+q =0\), then let \(y_p=e^{ax}[W_k(x)\cos β x+V_k(x)\sin β x]\), where \(k=\max(m, n)\).
(2) If \(k_1\) or \(k_2\) is equal to \(a+β i\), then let \(y_p=xe^{ax}[W_k(x)\cos β x+V_k(x)\sin β x]\), where \(k=\max(m, n)\).

Examples
(1) \(y''-y'=3e^{2x}+2x+1\). A particular solution of \(y''-y'=3e^{2x}\) is \(y_1=\frac{3}{2}e^{2x}\), and a particular solution to \(y''-y'=2x+1\) is \(y_2=-x^2-3x\). Thus, a particular solution to the equation is \(y_p=y_1+y_2=\frac{3}{2}e^{2x}-x^2-3x\), and the general solution is \(y=C_1e^x+C_2+\frac{3}{2}e^{2x}-x^2-3x\). Check ode;y_2-y_1-3*exp(2*x)-2*x-1;x ||.

(2) Solve \(y''-2y'+y=25e^{-x}\cos x\). The characteristic equation is \(k^2-2k+1=0⇒k=1\). Let a particular solution be \(y_p=e^{-x}(A\cos x+B\sin x)\). Substituting \(y_p'\) and \(y_p''\) in the DE and matching coefficients, we have \(A=3,B=-4\), \(y_p=e^{-x}(3\cos x-4\sin x)\), and \(y=e^x(C_1+C_2x)+e^{-x}(3\cos x-4\sin x)\) is the general solution to the DE. Check ode;y_2-2*y_1+y-25*exp(-x)*cos(x);x ||.

Examples ode;y_2+2*y_1-3*y-1;x || ode;y_2+4*y-x**2;x || ode;x;x_2+2*x_1+x-t*exp(t);t || ode;y_2+y_1-1-t-exp(t);t || ode;x_2+2*x_1+x-t/exp(t);t || ode;y_2+2*y_1+2*y-exp(x)*cos(x);x || ode;y_2-6*y_1+9*y-exp(-3*x);x || ode;y_2-4*y_1+5*y-x*exp(2*x);x ||.

Numerical Methods for Solutions of First-Order ODE

Slope fields Not all differential equations can be solved. In fact, many DEs are impossible to solve explicitly. But we can use a graphical approach to visualize their shapes of solution curves.

Consider a first-order DE of the form \(y'=f(x,y)\), where the function \(f(x,y)\) of \(x\) and \(y\) gives an explicit instruction on how to calculate the slope of the tangent line to the solution curve \(y\) at each point \((x,y)\). Or the curve of a solution to the DE has a tangent line at point \((x,y)\) whose slope is equal to the value of \(f(x,y)\).

If we draw a short segment of the tangent line with slope \(f(x,y)\) through each point \((x,y)\), then these line segments form a slop field or a direction field for the DE \(y'=f(x,y)\). Thus, each line segment gives a direction or slope for a solution curve passing through that point, and the patterns of the line segments display the shapes of all solution curves of the DE. Given initial conditions, we can connect the point \((x_0,y_0)\) of initial values with the other points following the directions of their tangent lines to form a particular solution curve.

Euler's method uses linear approximation to estimate a particular solution to a first-order DE. Suppose \(y=f(x)\) is a solution to a DE of the form \(y'=F(x,y)\) with an initial condition \(y_0=f(x_0)\). Then by linearization, the values of \(y=f(x)\) for \(x\) near the point \((x_0,y_0)\) can be approximated by \(f(x_0)+f'(x_0)(x-x_0)=y_0+F(x_0,y_0)dx\), in particular if \(y\) does not have any explicit form. The general procedure for Euler's method is as follows.

Starting from the point \((x_0,y_0)\), estimate \(y_1=f(x_1)≈ y_0+F(x_0,y_0)(x_1-x_0)\) for a new point \(x_1=x_0+dx\) near \(x_0\). Once \((x_1,y_1)\) is known, estimate \(y_2≈ y_1+F(x_1,y_1)(x_2-x_1)\) for \(x_2=x_1+dx=x_0+2dx\), where \(dx=\Delta x\). Continuing this process, we get \(y_n≈ y_{n-1}+F(x_{n-1},y_{n-1})(x_n-x_{n-1})\) for \(x_n=x_0+ndx\).

Example If \(y'=2x+y\) and \(y(0)=1\), estimate \(y(1)\). Set step size \(dx=0.1\). Then the nth iteration is given by \(x_n=0.1n, y_n=y_{n-1}+ 0.1(2x_{n-1}+y_{n-1})\). Starting from (0, 1), \(x_1=0.1,y_1=y(0.1)=1+0.1[2(0.1)+1]=1.12,y_2=y(0.2)=1.12+ 0.1[2(.2)+1.12]\) = 1.272. Continue this process, we get \(y(0.3)=1.4562,\cdots\), and \(y(1)≈4.1\). The exact solution is \(y=3e^x-2x-2\) and \(y(1)=3e-4\) ≈ 4.1548155. Check ode;y_1-2*x-y;x ||.

Errors of Euler's method Euler's method is not effective, but it provides a simple procedure to calculate a particular solution to a DE numerically. The errors of Euler's method are accumulating step by step. So the more steps it involves, the more errors are accumulated. To improve the accuracy, we can use smaller step size \(dx\). But keep in mind that using a very small step size may add rounding errors in approximating a solution, and how to reduce rounding errors is a topic in numerical analysis.

Motion with resistance Assume an object of mass \(m\) moves along a coordinate line with a position function \(s(t)\) and velocity \(v(t)\) at time \(t\). If the resisting force opposing the motion is proportional to velocity, then according to Newton's second law of motion, \(F=m\frac{dv}{dt}=-kv\) or \(\frac{dv}{dt}=-\frac{k}{m}v\) for \(k>0\). Solving the DE, we can find \(v(t)\) and \(s(t)\).

Examples This first order DE \(v'=-\frac{k}{m}v\) is separable, and it is similar to the exponential decay model. The general solution to the equation is \(v(t)=Ce^{-\frac{k}{m}t}\). Check ode;v_1+k*v/m;t ||. With initial condition \(v(0)=v_0\), we have \(v(t)=v_0e^{-\frac{k}{m}t}\). To find the position function \(s(t)\), solve the equation \(v(t)=\frac{ds}{dt}=v_0e^{-\frac{k}{m}t}\) for \(s(t)\) by integrating both sides with respective to \(t\). Thus, \(s(t)=\int v_0e^{-\frac{k}{m}t}=-\frac{mv_0}{k}e^{-\frac{k}{m}t}+C\) by ode;s_1-C*exp(-k*t/m);t ||. With initial condition \(s(0)=0,s(t)=\frac{mv_0}{k}(1-e^{-\frac{k}{m}t})\). As \(t\) tends to infinity, \(s(t)\) approaches \(\frac{mv_0}{k}\).

Examples In case of free-fall motion with air resistance taken into account, assume the force by air resistance is proportional to the velocity and acts opposite to the direction of falling. Choose the upward direction as positive direction. The net force \(F_g+F_r\) acting on an object is the sum of its weight \(F_g\) due to gravity and the force due to air resistance \(F_r\). According to Newton's second law of motion, \(F_g+F_r=-mg-kv=m\frac{dv}{dt}\), where \(k\) is a positive constant, or \(v'=-\frac{k}{m}v-g\). Solving this first-order linear DE, we have \(v(t)=-\frac{mg}{k}+Ce^{-\frac{k}{m}t}\), where \(C\) is any constant. Check ode;v_1+g+k*v/m;t ||.

As \(t\) tends to infinity, \(v(t)\) approaches \(-\frac{mg}{k}\), which is called a limiting terminal velocity. The terminal velocity can be simply obtained by letting the \(v'=0\), which means the acceleration is 0 or the velocity remains constant. In this case, \(-mg-kv=0\) implies \(v(t)=-\frac{mg}{k}\).

Logistic growth model A more realistic growth model considers constrains and effects of limited environmental resources (food and space) on the growth of a population. The logistic differential equation \(\frac{dy}{dt}=ky(1-\frac{y}{L})\) describes how growth rate is affected by the maximum population \(L\) (carrying capacity) for which the environment can sustain. It is similar to exponential growth model \(y'=ky\), but adds a factor \((1-\frac{y}{L})\) to adjust the growth rate. When \(y\) is small compared to \(L, 1-\frac{y}{L}≈ 1\) and \(\frac{dy}{dt}≈ky\), which means the population increases exponentially, in particular when the size of the population is relatively small.

If \(0< y< L\), the growth rate is positive and the population increases. But when \(y > L\), which means the population is greater than the maximum population the environment can support, the population starts to decrease, and the growth rate is negative. Both \(y=0\) and \(y=L\) are solutions to the equation because the growth rate is 0 for these two values. These two constant solutions are called equilibrium solutions.

By separation of variables, \(\int \frac{dy}{y(1-y/L)}=\int kdt=kt+C\). Note that the fraction \(\frac{1}{y(1-y/L)}=\frac{L}{y(L-y)}=\frac{1}{y}+\frac{1}{L-y}\). Substituting the fraction into the integration \(\int \frac{dy}{y(1-y/L)}=\int[\frac{1}{y}+\frac{1}{L-y}]dy=\ln|y|-\ln|L-y|=\ln|\frac{y}{L-y}|\)\(=kt+C\), or \(\ln|\frac{y-L}{y}|=-kt-C\), which implies \(|\frac{L-y}{y}|=e^{-C}e^{-kt}⇒ \frac{L-y}{y}=\pm e^{-C}e^{-kt}=Ae^{-kt}\), where \(A=\pm e^{-C}\). Therefore, \(y(t)=\frac{L}{1+Ae^{-kt}}\) and \(\displaystyle\lim_{t\to\infty}y(t)=L\).

Examples Verify the above results by ode;y_1-k*y*(1-y/L);t || ob by ode(y_1-k*y*(1-y/L),t) ||.

Orthogonal trajectory An orthogonal trajectory of a family of curves is a curve that intersects each member curve of the family at right angles.

For example, each member of the family \(y=mx\) of straight lines is an orthogonal trajectory of the curves of the family \(x^2+y^2=r^2\). A line passes the origin and intersects each circle at a right angle. The two families are orthogonal of each other.

To find the orthogonal trajectories of a family of curves defined by the equation \(y'=f(x,y)\), we need to solve the DE \(y'=-\frac{1}{f(x,y)}\), because the product of the slopes at any intersection point \((x,y)\) of curves from the two families is -1.

Examples
(1) To find the orthogonal trajectories of the straight lines \(y=mx\) or \(y'=m=\frac{y}{x}\), we need to solve the equation \(y'=-\frac{x}{y}\). By separation of variables, \(\int yy'dx=\int -xdx⇒\frac{y^2}{2}=-\frac{x^2}{2}+C\) or \(x^2+y^2=r^2\) for \(r^2=2C\). Check ode;y_1+x/y;x ||.

(2) Show the orthogonal trajectory of parabolas \(y=ax^2\) is a family of ellipses. Since \(y=ax^2,y'=2ax,a=\frac{y}{a^2}\). The orthogonal trajectory satisfies \(y'=-\frac{1}{2ax}=-\frac{x}{2y}⇒2yy'=-x⇒∫2ydy=-∫xdx⇒y^2=-\frac{x^2}{2}+C\)\(⇒2y^2+x^2=2C\) by ode;y_1+x/(2*y);x ||.

Newton's law of cooling/warming Let \(y(t)\) be the temperature of an object at time \(t\), and \(T_e\) be the temperature of the environment. Then Newton's law of cooling can be described by the DE \(\frac{dy}{dt}=-k(y-T_e)\), where \(k>0\) is the constant of proportionality. Thus, if \(y > T_e, y'< 0\), or the rate of change in temperature is negative, which means the object's temperature is decreasing.

Similarly, if \(y < T_e, y'>0\), which means the object's temperature is increasing. By separation of variables, \(\int \frac{dy}{y-T_e}=-\int kdt=-kt+C⇒\ln|y-T_e|=-kt+C⇒ |y-T_e|=e^Ce^{-kt}⇒ y=T_e+Ae^{-kt}\).

Example An object at 80°C is placed at a room where temperature is 25°C. After 4 minutes later, the temperature has dropped to 60°C. Find the time when the object is at 50°C. Let \(y(t)\) be the temperature of the object at time \(t\). Since \(y'=-k(y-T_e),T_e=25,y(0)=80⇒55=Ae^{-k(0)}⇒A=55,y(4)=60=25+55e^{-4k}\) ⇒ \(k\) ≈ 0.113. Thus, \(y=25+55e^{-0.113t}\) and \(50=25+55e^{-0.113t}⇒t≈6.98\) minutes. Check the results by ode;y_1+k*(y-80);x || slv;35-55.0*exp(-4.0*k);k ||.

Mixture problem Suppose a chemical in a liquid solution flows into a container that is already filled with a solution of a known amount of the same chemical dissolved into it. During the time the chemical is added to the container at a given rate, the liquid is thoroughly mixed in the container by stirring and is also drained from the container at a known rate. As time goes on, the amount of chemical in the container changes, and in general the rate of change is \(\frac{dy}{dt}\) = (rate in) - (rate out). The question of interest is how to determine the amount of the chemical at a given time.

Assume at time \(t\) the amount of the chemical is \(y(t)\), and the volume of the liquid in the container is \(V(t)\). Then the concentration at time \(t\) is \(\frac{y(t)}{V(t)}\). If the outflow rate is \(R_o\) volume per unit time, then the rate out is \(\frac{y(t)}{V(t)}R_o\). Meanwhile, the rate in is \(CR_i\), where \(C\) is the concentration of the inflow liquid, and \(R_i\) is the rate of the liquid entering the container. Now the rate of change of the chemical is \(CR_i- \frac{y(t)}{V(t)}R_o\) at time \(t\). Use a differential equation \(\frac{dy}{dt}=CR_i-\frac{y}{V(t)}R_o\) to describe this process, solve the equation, and we can obtain the amount of the chemical at a given time \(t\).

Example A tank contains 64 gallon of a mixture in which 16 lb of salt are dissolved. A solution containing 2 lb salt per gallon flows into the tank at a rate of 4 gal/min. The solution is thoroughly mixed by stirring, and drains from the tank at a rate of 1 gal/min. Determine the amount of salt at time \(t\) = 20 minutes. Let \(y(t)\) be the amount of salt at \(t\). Then the rate of change of salt is \(\frac{dy}{dt}\). The rate of salt entering in the tank is \(CR_i\) = 2(4) = 8 lb/min. The concentration of salt is \(\frac{y(t)}{64+(4-1)t}\) lb/gal, the rate of the liquid exits the tank is \(R_o=1\) gal/min, so the rate of salt exits the tank is \(\frac{y}{64+3t}R_o=\frac{y}{64+3t}\) lb/min, and the differential equation is \(\frac{dy}{dt}=8-\frac{y}{64+3t}\). The general solution of the equation is \(y=C(64+3t)^{-\frac{1}{3}}+6t+128\) by ode;y_1-8+y/(64+3*t);t ||. Given \(y(0)\) = 16 ⇒ \(C\) = -448, and \(y(20)\) ≈ 158.16 lb.

Partial differential equation: [pde; expr; iv1; iv2], [pde(expr, iv1, iv2)]

A first-order partial differential equation (PDE) is an equation that involves an unknown function of \(n\) variables and its first partial derivatives. The first order PDE takes the form \(f(x_1,x_2,\cdots,x_n,y'_{x_1},y'_{x_2},\cdots,y'_{x_n})=0\). The unknown function \(f\) is multivariate and has at least two independent variables. If all coefficients of the unknown function and its partial derivatives are constants or functions of the independent variables only, the PDE is also linear.

For example, if \(n=2\), the unknown function \(z=f(x,y)\) in the PDE is a bivariate function, the PDE such as \(x\frac{∂z}{∂y}-y\frac{∂z}{∂x}=0\) is a first-order linear PDE, and the general solution is in form of \(z=F(x^2+y^2)\). The equation \(xz_x+yz_y=z\) is first-order linear PDE, and its general solution has the form \(z=F(\frac{y}{x})x\). In these examples, the unknown function \(z\) is bivariate and involves just two independent variables \(x\) and \(y\). In this section, we only focus on solving first order, linear PDE that involves only two independent variables by the program.

Use the code pattern "pde; expr; iv1; iv2" to find the solution to some form of the first-order linear PDE of two variables, where "pde" (partial differential equation) is the operation name, "expr" (without the portion "=0") is the expression of the partial DE, and "iv1; iv2" the independent variables involved in the PDE.

A valid expression "expr" of a PDE must be the first order partial derivatives, and a valid partial derivative such \(\frac{∂z}{∂x}\), the partial derivative of \(z(x,y)\) with respect to \(x\), must be written as "z_x". Similarly "z_y" represents \(\frac{∂z}{∂y}\), the partial derivative of \(z\) to \(y\). Write the name of an unknown function, followed by an underscore and the name of either independent variable. Write first order partial derivatives using other variable names in a similar fashion.

Attention Use the code pattern "pde(expr, iv1, iv2)" to get the same results as "pde; expr; iv1; iv2".

Attention Use 'dif" module to write the a partial derivatives such as "dif(f(x,y),x)" for f_x, and "dif(f(x,y),y)" for f_y. Then apply "pde" module to get the general solution for f(x, y). Combining "pde" and "dif" modules for solving a partial DE is conceptually straightforward. You must write the unknown function as f(x, y) and write variable x or y as a single letter. Don't mix using derivatives such as "f_x" and "dif(f(x, y),x)" for f_x in the same DE, or "f_y" and "dif(f(x, y),y)" for f_x. The partial DE written in this way clearly specifies each term of partial derivatives, unknown functions, and variables, making it ready for using the "pde" module.

Examples
(1) Type pde;x*z_x+y*z_y-z;x;y ||, and it gives the general solution \(z=F(x^2+y^2)\) to the PDE \(x\frac{∂z}{∂x}+y\frac{∂z}{∂y}=z\), where \(z\) represents the unknown function \(z(x,y)\), and "x" and "y" the two independent variables.
(2) The code pde;u_x-u_y;x;y || gives the general solution \(u=F(-x-y)\) to the PDE \(\frac{∂u}{∂x}-\frac{∂u}{∂y}=0\).
(3) If f_x = 0 for all (x, y), then f(x, y) is by pde;f_x;x;y || or pde(f_x,x,y) || or pde(dif(f(x,y),x)) ||.
(4) If f_x = g(x) for all (x, y), use pde;f_x-g(x);x;y || pde(dif(f(x,y),x)-g(x)) || for f(x, y). If f_x = g(y), use pde;f_x-g(y);x;y || pde(dif(f(x,y),x)-g(y)) || for f(x,y).

Example pde;z_x+y;x;y || pde;z_y-x;x;y || pde;f_x-x*y;x;y || pde;z_x-x**2*y;x;y || pde;z_y-x*y**2;x;y || pde;z_y-z_x-1;x;y || pde;z_x-z*x*y;x;y || pde;x*z_y-y*z_x;x;y || pde;w_u-w_v;u;v || pde;w_x+w_y;x;y || pde;y*z_x-y-z;x;y || pde;x*z_y-z;x;y || pde;z_x-z_y-3*z;x;y || pde;x**2*z_x+y**0.5*z_y-1;x;y || pde;x*y*z_x+y**2*z_y-x*y;x;y || pde;-2*z_x+4*z_y+5*z-exp(x+3*y);x;y || pde;z_x-2*x*z_y;x;y || pde;y*z_x-x;x;y || pde;y*z_x-z;x;y || pde(y*z_x-z,x,y) || pde(y*z_x-x,x,y) ||.

Ordinary first-order linear differential equation systems: [ods; iv; eq1; eq2; ... ], [ods(iv,eq1, eq2,...)]

The growth rate of one population (prey) may affect the rate of another population (predator) that lives in the same environment. This situation leads to two differential equations that must be solved simultaneously in order to determine the two populations. In this section, we only focus on solving first order linear system of ODEs.

In general, the first-order system of ODEs takes the form \(y_i'=f_i(x,y_1,y_2,\cdots,y_n)\) for \(i=1,2,\cdots,n\) for the independent variable \(x\) on some interval. The right side of each equation in the system contains no derivatives. Just as in a single ODE, if the coefficients of the unknown functions and its derivatives of each equation in the system are constants or functions of \(x\) only, then the system is first-order linear system of ODEs.

For example, \(y_1'=2xy_1+y_2\), \(y_2'=3x-2y_1\) is a first-order linear system of ODEs. A higher order ODE or ODE system may be made into a first-order system by change of variables.

Examples
(1) \(y''-2y'+5y=0\) can be written as a first-order system by letting \(t=x,u(t)=y(t),v(t)=y'(t)\). Then \(u'=v,v'=2v-5u\). Check ode;y_2-2*y_1+5*y;x || ods;t;u_1-v;v_1-2*v+5*u ||.

(2) Write \(y'''-3y''-2y'+4ty=0\) as a system by letting \(u(t)=y,v(t)=y',w(t)=y''\). Then \(u'=v,v'=w\),\(w'=3w+2v-4u\). Check ods;t;u_1-v;v_1-w;w_1-3*w-2*v+4*u ||.

Use "ods; iv; equ1; equ2; ... " to help solve a system of first-order linear ODEs, where "ods" (ordinary equation system) is the operation name, "iv" the independent variable involved in the system, and "equ1; equ2; ..." represents each ODE (without the portion "= 0") in the system.

Attention Also use function call "ods(iv,eq1, eq2,...)" for solutions to first-order linear ODEs, and you get the same results as "ods; iv; equ1; equ2; ... ".

Examples
(1) Type ods;t;x_1-3*x+2*y;y_1-2*x+y ||, and you get the solution to the system \(x'=3x-2y, y'=2x-y\), where \(t\) is the independent variable (dummy) "iv", and \(x(t)\) and \(y(t)\) are the unknown functions. The general solution is \(x(t)=2Ate^t+(A+2B)e^t,y(t)=2Ate^t+2Be^t\). If \(x(0)=1;y(1)=4e\), then \(A=3,B=-1\), and the particular solution is \(x(t)=(6t+1)e^t,y(t)=(6t-2)e^t\). Also by ods(t,x_1-3*x+2*y,y_1-2*x+y) ||.

(2) ods;t;x_1-y-z;y_1-x+z;z_1-x-y || gives the solution to the system \(x'=y+z,y'=x-z,z'=x+y\).

Examples ods;t;x_1-x-2*y;y_1-2*x-3*y || ods;t;x_1+x-y;y_1-x+y || ods;t;x_1-2*y;y_1-3*x || ods;t;x_1-y_1+2*x-3*y;y_1-2*x_1+x-2*y || ods;t;x_1-2*x-3*y+w;y_1-x+y;z_1-y+w;w_1-2*x || ods;t;x_1+x+y-exp(-t);y_1-x-y-exp(-t) ||.