DSolve
 Introduction

Let us examine the Mathematica function DSolve which is used to find analytical solutions of first order ordinary differential equations.  

The general form of the DSolve command is

                         DSolve[equation, dependent variable, independent variable]

We will first learn how to specify an equation for DSolve.  Then we will try it out, thereby learning something
about the form of the output from running DSolve (which is a replacement rule).  After showing how to get a function rather than a replacement rule as an output from DSolve, we carry out a number of examples.  We conclude with examples of how to automate the production of tables or a graph sequence from solutions of the equation.

Specifying an Equation for DSolve

The syntax required by DSolve for a differential equation is very similar to standard mathematical notation.  Derivatives are denoted by a prime, which is typed as the apostrophe character.  Multiple derivatives
are denoted by repeated apostrophes.  

Dependent variables always exhibit explicitly the independent variable as a functional argument.  

Equality in the equation is denoted by the double equals sign ( ==).  Here is an example, typed
first in standard mathematical notation:

(dy/dx) +  y = Sin(3x)        where y is a function of x

and now  in a form suitable for DSolve,

					y'[x] +  y[x] == Sin[3x]  
					
The second argument of DSolve, the dependent variable (y), is written with the functional argument as y[x],

and the independent variable as x.  Define this DE as an equation   deq = ........

Using what know about Solve   or look up DSolve and get DSolve expression ready to go.

ClearAll["`*"]
deq = y'[x] + y[x]==Sin[3 x]
sol1 = DSolve[deq,y[x],x]

Without a (ONE) BC we get a constant of integration.

Extract your solution

y[x_] = y[x]/.sol1[[1]]/.C[1]-> cc1  (* need to do this C[1] is protected  cc1 IS NOT! *)

After the fact:  introduce a BC to remove C[1];   Let y[0]==1 for an example

Find y[x] “the hard way” and plot it around x = 0 (where you know the BC)

Normally we would start with sol2 = Solve[y[0] ==1,cc1], do the extraction and if necessary N[...]  and possibly Simplify[....]

cc1 = cc1/.Solve[y[0] ==1,cc1][[1]]//N
y[x]

Plot[y[x],{x,-1,1}]

Now merge all the input cells, find and plot y[x] introducing the BC into DSolve (just like we did with Solve).  TRIGGER ALERT:  You’ll have to do some editing

ClearAll["`*"]
deq = y'[x] + y[x]==Sin[3 x];
y[x_] = y[x]/. DSolve[{deq, y[0] == 1},y[x],x][[1]]
(*  TEST -- is y[x] a solution??? *)

y'[x] + y[x]==Sin[3 x]
Plot[y[x],{x,-1,1}]

Vertical motion with a v dependent drag force (~ (v^2))  

first, find velocity  v[t] assuming vo is up

acceleration a  = v’[t] = (F/m) (-g + k ((v[t])^2))

also find to = time to reach v = 0  (at top of motion).

Going UP:

ClearAll["`*"]
deq = v'[t] == -g - (k/m)(((v[t]))^2);
sol = DSolve[{deq,v[0]==vo},v[t],t];
v[t_] = v[t]/.sol[[1]]
to = t/.Solve[v[t]==0,t][[1]]

What does ArcTanh look like      to ~ ArcTanh[CC vo]     Let CC vo  = z

Plot[ArcTan[z],{z,0,5}]

Chose some values and plot;  we will chose m = 1; k = 1; g = 1;   THEN (Sqrt[(k/(m g))])  = 1 
NOTE that to = ArcTan[vo]

m = 1; k = 1; g = 1; vo =1.5;
to
v[t]
vplot = Plot[v[t],{t, 0,to},PlotLabel->Style[Framed["v[t] and (Subscript[t, o]) (red dot) TOP OF MOTION"],16,Blue,Background->Lighter[Yellow]]];
topoint =Graphics[{PointSize[Medium],Red,Point[{to,0}]}];
Show[vplot,topoint]

Now find x[t]:    a = x’’[t]       v[t] = x’[t]

Evaluation-> QuitKernal-> Local;
ClearAll["`*"]
deq = x''[t] == -g - (k/m)(((x'[t]))^2);
sol = DSolve[{deq,x[0]==0,x'[0] == vo},x[t],t]
x[t_] = x[t]/.sol[[3]]
v[t_] = x'[t];
to = t/.Solve[{v[t]==0},t][[1]]

m = 1; k = 1; g = 1; vo =1.5;
xplot = Plot[x[t],{t, 0, to},PlotLabel->Style[Framed["x[t] and (Subscript[t, o]) (red dot) where v = 0"],16,Blue,Background->Lighter[Yellow]]];
topoint =Graphics[{PointSize[Medium],Red,Point[{to,x[to]}]}];
Show[xplot,topoint]
to

Log[-2]

Some more examples:

y'[x]+y[x]==1   no BCs

ClearAll["`*"]
sol = DSolve[y'[x] + y[x] == 1, y[x], x]
y[x_] = y[x]/.sol[[1]]

Here's an "exact" equation:      x y’[x] + 3x + y[x] == 0  (single solution without any BCs)

First by hand:
rewrite in M N form:     x dy + (3 x + y) dy = 0

(∂_y)u = x;        (∂_x)u = (3x + y)      
Integrating  u = x y + g[x]           u = 3 ((x^2)/( 2)) + x y +h[y]
SAME if g[x] =  3 ((x^2)/( 2))  and h[y] = 0
u[x,y] =  3 ((x^2)/( 2)) + x y   cc.
Solve for y[x] = (cc/x) - (3/2)x   
NOW use DSolve:

Remove["`*"]
sol = DSolve[x y'[x] + 3x + y[x] == 0, y[x], x]/.C[1]-> cc
y[x_] = y[x]/.sol[[1,1]]

Here's another example.  It isn't "exact", but because the left hand side is linear in y it could be solved by finding an integrating factor depending only on x  (I have no idea if Mathematica uses that approach).

x y’[x] + (1 + x)*y[x] == Exp[x]

Remove["`*"]
y[x_] = y[x]/.(DSolve[x y'[x] + (1 + x)*y[x] == Exp[x], y[x], x]/.C[1]-> cc)[[1,1]]

This DE gives us two integration constants because of the y’’:

x^2 y''[x]+x y'[x]+x^2 y[x]==0

Remove["`*"]
DSolve[x^2 y''[x] + x y'[x] + x^2 y[x] == 0, y[x], x]

Try the above with two BCs that you think will yield a solution; Try plotting it.

Remove["`*"]
y[x_] = y[x]/.(DSolve[{x^2 y''[x] + x y'[x] + x^2 y[x] == 0, y[0] == yzero, y'[0] == 0   }, y[x], x])[[1]]

yzero = 10;
Plot[y[x],{x,0,20}]

Here's another example where we get two solutions but single integration constant
y’[x]^2 + y’[x] - 2  == 0

Remove["`*"]
DSolve[((y'[x])^2) + y'[x] - 2  == 0, y[x],x]

Here's another exact equation with a boundary condition specified.  Notice that the boundary condition must be used "==", not "=".

DSolve[{x y'[x]/y[x] + Log[y[x]] == 2, y[1] == 1},
       y[x], x]

Note that the solution is complete, without any arbitrary constants remaining:

y[x]/.%[[1]]

If we add a BC at y[0]  M chokes  --  (y would blow up)

DSolve[{x y'[x]^2 +  y'[x] - x y[x]  == 0,y[0] ==0}, y[x],x]

Normally, if M cannot find the solution, it just gives back your input.

RL Circuit
                                                 
Consider the above circuit where a DC Voltage, V = Vo has been applied for a long time and at t = 0 is switched to a short (V = 0).

(Subscript[I, o]) =  ((Subscript[V, o])/R) and for t>0, Kirchoff’s Law (V around loop) gives 
                                       L(dI/dt) + I R = 0;

ClearAll["Global`*"]
Clear[Vzero]
sols =DSolve[{L*II'[t]+ R* II[t]==0,II[0] == (Vzero/R)},II[t],t] 
y[t_] = II[t]/.sols[[1]]; (* This trick extracts the solution; NOTE = rather than := in definition *)
Print["Here is I[t]: I[t] = ",Simplify[ y[t]]]

Clearly, this is an exponential decay from Izero = ((Subscript[V, o])/R) to zero; if R is increased, the time constant increases and if L increases the time constant decreases.  Here is a Manipulate plot with R and L varying independently.

Vzero = 10;
Manipulate[{{R, L},Plot[(((E^(-((R t)/L))) Vzero)/R),{t,0,50},PlotRange->{{0,50},{0,100000}},AxesOrigin->{0,0}]},{R,0.0001,.005,0.00001},{L,0.0001,1,0.0001}]

Note that when you increase R, the amplitude changes (as (Vzero/R)). 
Now we look at the total energy delivered to the resistor:   P = (I^2)R

Print[Style["  Here is the voltage
 i) across the Resistor = I*R AND 
ii) across the capacitor = Q/C : ",16,Blue, Bold]]
Print["    ",Simplify[y[t]]*R,"             ",Simplify[ y[t]]/C]

TextCell["Below we plot the current vs. time for fixed Vzero and R; we vary C (which varies τ = RC (the time constant for the decay:","Subsection"]
Vzero = 10; R = 2.5;
Manipulate[{C,Plot[(((E^(-(t/(C R)))) Vzero)/R),{t,0,30},PlotRange->All,AxesOrigin->{0,0}]},{C,0.001,5,0.01}]
TextCell["Note the as you increase C, which increases τ = RC (the time constant for the decay), the current falls slower and slower.","Subsection"]
Print[]
TextCell["Now we plot the Voltage drop across the resistor == I*R, and the voltage across the capacitor == Q/C; THEY ARE EQUAL -- look at the circuit (I've displaced I*R slightly downwards so you can see it).  Again,  τ = RC (the time constant appearing in both functions), the changes in V_R and V_C are slower and slower for larger C (larger τ = RC).","Subsection"]
Manipulate[{C,Plot[{0.9(E^(-(t/(C R)))) Vzero,(E^(-(t/(C R)))) Vzero},{t,0,30},PlotRange->All,AxesOrigin->{0,0}]},{C,0.001,5,0.01}]
TextCell["Remember: these are the same -- one is displaced so you can see it!","Subsection"]

The general solution to the Series RLC circuit.  

First we look at under, critical,and over damped separately.

Here is Underdamped:

ClearAll["Global`*"]
Print["This is underdamped ONLY -- meaning that β < ωo"]
Imax = 1;
ω1 = (Sqrt[((ωo)^2) - (β^2)])
Manipulate[Plot[Imax *(E^(-β t)) (Cos[ (Sqrt[((ωo)^2) - (β^2)]) t] + (β/(Sqrt[((ωo)^2) - (β^2)]))Sin[ (Sqrt[((ωo)^2) - (β^2)]) t]),{t,0.1,30},PlotRange->1],{β,0,5,.1},{ωo,1,5}]

Let’s find critically damped solution by taking the limit of underdamped as ωo -> β, which is the same thing as letting ω1 ==(Sqrt[((ωo)^2) - (β^2)])  -> 0

ClearAll["Global`*"]
Limit[Imax*(E^(-β t)) (Cos[ ω1 t] + (β/(ω1))Sin[ ω1 t]),ω1->0]

ClearAll["Global`*"]
Print["This is critially damped ωo = β "]
Imax = 1;
Manipulate[Plot[(E^(-t β)) Imax (1+t β),{t,0.1,30},PlotRange->1],{β,0.001,5}]

i(t) = Izero{Exp[-βt] * {cosh (Subscript[ω, 2])t  + (β/(Subscript[ω, 2])) sinh (Subscript[ω, 2])t }
(Subscript[ω, 2]) = (Sqrt[(β^2)-(Subscript[ω, o]^2)])

ClearAll["Global`*"]
Print["To Show what Cosh and Sinh Look Like (Sinh passes through {0,0}:"]
Plot[{Cosh[x],Sinh[x]},{x, 0,2}]

Print["This is overdamped ONLY -- meaning that β > ωo"]
Imax = 1;
Manipulate[Plot[Imax *(E^(-β t)) (Cosh[ (Sqrt[ (β^2) -((ωo)^2)])  t] + (β/((Sqrt[ (β^2) -((ωo)^2)]) ))Sinh[ (Sqrt[ (β^2) -((ωo)^2)])  t]),{t,0.1,30},PlotRange->1],{β,0.20001,2},{ωo,.2,2}]

Do we get same critical damped solution from OVERDAMPED solution if we let ωo -> β?
This is what we got for critical damped I(t):   (e^(-t β)) Imax (1+t β)

Clear[Imax]
Limit[Imax *(E^(-β t)) (Cosh[ (Sqrt[ (β^2) -((ωo)^2)])  t] + (β/((Sqrt[ (β^2) -((ωo)^2)]) ))Sinh[ (Sqrt[ (β^2) -((ωo)^2)])  t]),ωo -> β]

Now we just let M do the whole job:

ClearAll["Global`*"]
sol = DSolve[{y''[t]+2 β y'[t]+  ωo^2  y[t]==0,y[0]==Imax,y'[0]==0},y,t]
β = 1/2; ωo = 1; Imax = 1;

Imax = 1; 
Manipulate[{Plot[(1/(2 (Sqrt[(β^2)-((ωo)^2)])))Imax (-(E^(t (-β-(Sqrt[(β^2)-((ωo)^2)])))) β+(E^(t (-β+(Sqrt[(β^2)-((ωo)^2)])))) β+(E^(t (-β-(Sqrt[(β^2)-((ωo)^2)])))) (Sqrt[(β^2)-((ωo)^2)])+(E^(t (-β+(Sqrt[(β^2)-((ωo)^2)])))) (Sqrt[(β^2)-((ωo)^2)])),{t,0.1,30},PlotRange->1]},{β,0,5},{ωo,1,5}]

Is the behavior for imaginary args for under and over damped solutions OK?  are these solutions equivalent?  

There is a connection between complex hyperbolic and complex trigonometric functions:

Cosh[z] == Cos[I*z]
Sinh[z] == -I * Cos[I*z]
Cos[z] == Cosh[I*z]
Sin[z] == (-I)*Sinh[I*z]
Sinh[z] == I*Sin[I*z]
Sin[z] == I*Sinh[I*z]

Cosh(z) = Cos(iz)

Sinh(z) = - iCos(iz)

Cos(z) = Cosh(iz)

Sin(z) = - iSinh(iz)

Sinh(iz) = iSin(z)

Sin(iz) = iSinh(z)

Let’s look at y[x] vs. y

ClearAll[x,y,yy];
eqn = {y'[x]==y[x]^2, y[0]==1}
sol= DSolve[eqn, y[x], x]
sol = Flatten[sol]
yy[x_] = y[x]/.sol[[1]]
yy[x]
Plot[yy[x],{x,0,1}]

eqn = {y'[x]==y[x]^2, y[0]==1};
sol= DSolve[eqn, y, x]
sol = Flatten[sol]
yy[x_] = y[x]/.sol[[1]]
yy[x]
Plot[yy[x],{x,0,1}]

Find solutions to the equation y’[x] == 3 y[x]

ClearAll["Global`*"]
DSolve[y'[x]==3 y[x],y[x],x]

stipulate an initial condition   I chose y[2]  = π

(*stipulating an initial condition*)
ClearAll["Global`*"]
DSolve[{y'[x]==3y[x],y[2]==π},y[x],x]
sol = Flatten[%]
y[x_] = y[x]/.sol[[1]]

Now plot y[x];  Use Mouse to confirm BC

Dynamic[MousePosition["Graphics"]]
Plot[y[x],{x,0,2}]

Confirm that y[2]

(*confirming that y[2]=π  Here we use the Replacement Operator twice i) like above; ii) x=2 *)
Print["(confirming that y[2]=π):"]
y[2]

Confirm that y[x] is solution to DE  using a (==) test

(*confirming that y'[x]=3 y[x] Again we use the Replacement Operator twice i) like above for LHS y[x]; ii) and again for RHS y[x] *)
Print["(confirming that y'[x]=3 y[x]):"]
y'[x] == 3 y[x]

If we do not specify initial conditions then the general solution will involve an arbitrary constant. Below, we use Table to create a list of particular solutions and then plot all of them at once.  Notice that this differential equation produces two different solutions so that the output from DSolve is a list with two elements.  The solutions are nearly the same: a plus sign is switched to a minus sign to go from the first to the second solution.  

Also, observe the construction y[x]/. sol/.C[1]→k in the Table function. We first replace y[x] with the solution and then replace C[1] with k before letting k run from −10 to 10 in steps of 1. 

It might seem that Table[y[x]/.sol, {C[1], −10, 10, 1}] would be simpler, but unfortunately this will not work. Mathematica complains that “Tag C in C[1] is Protected.” We also use Tooltip so that as the mouse is moved over different solutions the value of the constant k is displayed.

(*plotting different solutions*)
(*mouseover displays C[1]*)

sol=DSolve[y'[x]==(x^2-y[x])/(x+y[x]),y[x],x]
(*vary C[1] to produce list of solutions*)
sol2 = sol/.C[1]->k
particularSolutions=Table[y[x]/.sol2,{k,-10,10}];
(*plot all solutions in list*)
Plot[particularSolutions,{x,-4,4}]

Mathematica can solve higher order differential equations also.   The next example solves a fairly simple equation; first with no BCs and then with given y[0] and y’[0].

 Solve:   y’’[x] + y[x] == x;  no conditions.    Then solve again letting y[0] == 2, y’[0] == 0.  Plot this solution

ClearAll["Global`*"]
(*solving a second order differential equation*)
DSolve[y''[x]+y[x]==x,y[x],x]

Now solve it again letting y[0] == 2, y’[0] == 0.  Plot this solution

(*solving a second order differential equation with initial conditions for both y and y'*)sol = DSolve[{y''[x]+y[x]==x,y[0]==2,y'[0]==0},y[x],x]
Plot[y[x]/.sol,{x,-20,20}]

Consider a system of first order differential equations:
x’ [t]  == x[t] - 10y[t]
y’[t]  == 15x[t] + y[t]

We can do this with DSolve. We list both differential equations as the first argument to DSolve.  Notice that the second argument is also a list—a list of the two unknown functions.

ClearAll["Global`*"]
(*solving a pair of first order linear equations simultaneously*)DSolve[{x'[t]==x[t]-10y[t],y'[t]==15x[t]+y[t]},{x[t],y[t]},t]
sol = Flatten[%]
x[t_] = x[t]/.sol[[1]]
y[t_] = y[t]/.sol[[2]]
(*finding the initial values in terms of the constants C[1] and C[2]; turns out {x[0], y[0]}== {C[1], C[2]}*)
Print[]
Print["Initial values {x[0],y[0]} in terms of the constants C[1] and C[2]: "]
Print["          {x[0],y[0]} == ",{x[0],y[0]}]

(*using Locator to Manipulate the initial conditions while plotting the solution*)
Manipulate[ParametricPlot[
{x[t],y[t]}/.
{C[1]->pt[[1]],C[2]->pt[[2]]},
{t,0,length},PlotRange->{{-10,10},{-10,10}}   (* changing the tmax for the problem using controls below *)
],
(*manipulate controllers*)
{{length,.5,"Length"},0.001,4},
{{pt,{RandomReal[],RandomReal[]}},Locator}
]

DSolve can also solve differential-algebraic equations. The syntax is the same as for a system of ordinary differential equations.

Remove["`*"]
eqns = {f''[x] == g[x], f[x] + g[x] == 3 Sin[x], f[π]==1, f'[π] == 0};
sol = DSolve[eqns, {f[x],g[x]},x]
f[x_] = f[x]/.sol[[1]]
g[x_] = g[x]/.sol[[1]]

Let’s plot the solutions substitued into their sum to see if they satisfy the algebraic relation:  i.e, test f[x] + g[x] == 3 Sin[x]

Plot[{f[x],g[x], f[x] + g[x],3 Sin[x]+.5}, {x,-5,5}]

Consider a system of ODEs. 

Each solution is labeled according to the name of the function (here, x and y), making it easier to pick out individual functions.  (this takes a long time to evaluate):

ClearAll["Global`*"];
eqns={((t^2)+1)*x'[t]==-t*x[t]+y[t]-Sign[t], ((t^2)+1)*y'[t]==-x[t]-t*y[t]+t*UnitStep[t], x[0]==-1/2, y[0]==2};
sol=DSolve[eqns, {x[t],y[t]},t]

((t^2)+1)  x'[t] + t*x[t] - y[t]+Sign[t] = 0;      ((t^2)+1)*y'[t]+x[t]+t*y[t]-t*UnitStep[t] ==0

This plots the solutions.

Plot[{x[t]/.sol, y[t]/.sol}, {t,-10,10},PlotLegends->"Expressions"]

Test x[t] with the two solutions (one containing ArcTan,  the other t):

Plot[((-1+4 t+2 (ArcTan[t])+2  ((1/2) Log[1+(t^2)]))/(2 (1 +(t^2)))),{t,-5,0}]

Plot[((-1+4 t+2 (-t)+2 t (0))/(2 (1 +(t^2)))),{t,-5,0}]

They are the same;  assume that the two solutions for y are also the same.

({{DSolve[eqn,u[x,y],{x,y}], solve a partial differential equation for u[x,y]}})
Write with pretty formating the diffyque below here:

(x^2)  ((∂ u(x,y))/(∂ x)) + (y^2)  ((∂ u(x,y))/(∂ y))  - (x + y) u(x,y) = 0

Convert to Traditional Form:

(x^2)  ((∂ u(x,y))/(∂ x))+(y^2)  ((∂ u(x,y))/(∂ y))-(x+y) u(x,y)=0

Clear[x,y];
eqn=x^2*D[u[x,y],x]+y^2*D[u[x,y],y]-(x+y)*u[x,y];
sol=DSolve[eqn==0, u, {x,y}]

Here is the general solution to a linear first-order PDE. 

In the solution, C[1] labels an arbitrary function of ((-x+y)/(x y)).

fn=u[x,y]/.sol[[1]]
fn/.{C[1][t_]-> Sin[t^2] +(t/10)}
u[x_,y_] = fn

Heat Diffusion Equation --One Dimensional

Consider a rod with some initial temperature profile  as a function of it’s length; the behavior with time DE is given by:

deqn=D[u[x,t],t] == d D[u[x,t],x,x];
TraditionalForm[deqn]

For a initial T profile we put in a rectangular distribution of width w at the center of the rod

g[x_,w_]:= Piecewise[{{0,x< -w},{100,-w ⩽x⩽+w},{0,x>w}}]
w = .05;
Plot[g[x,w],{x,-1,1}]

d=0.1;
sol=NDSolve[{deqn,u[x,0] == g[x,1],u[-10,t] == 0,u[10,t] == 0},u,{x,-10,10},{t,0,100}]

Plot3D[Evaluate[u[x,t]/.sol[[1]]],{x,-10,10},{t,0,100},PlotPoints ->40,Mesh ->False,PlotRange ->All]

There are three different forms:

({{DSolve[eqn, y, x], ODE}, {DSolve[{eqn1,...,eqnM}, {y1,...,yN}, x], simultaneous ODE}, {DSolve[eqn, y, {x1, x2, ...}], PDE}})

The first form, DSolve[eqn, y, x], is the normal usage for ODEs.  This description writes the second argument as y, but in fact you can use the y[x] form, i.e. DSolve[eqn, y[x], x] .  If you use the plain "y" form you get a "pure" function for the result, as we'll explore later on.

Also we can specify more than one equation,  DSolve[{eqn1, eqn2}, y[x], x], allows you to give some boundary conditions, as we'll see.

The second form, DSolve[{eqn1, ..., eqnM}, {y1, ..., yN}, x] shows that you can specify simultaneous ODEs (or a list of ODEs).  Specifying N more than one dependent variables and M equations (with M≥N) lets you solve a set of N simultaneous ODEs with the same independent variable x.

The third form,  DSolve[eqn, y, {x1, x2, ...}], shows that you specify  more than one independent variable, which would let us solve a PDE.