Series Solution Examples

This notebook is by Steven Amgott.  Please send any questions or comments to samgott1@swarthmore.edu.  Feel free to use and distribute this notebook, but keep this author information in any copy you use or distribute.

In any input cell containing xxx , you must replace it with your input before evaluating the cell.  In general, anything in magenta is something you can, and possibly should, change (in particular if you wish to do a different example).

Initialization  (Can be skipped, if you answer "Yes" to the initialization request.)

The cell in this section is an initialization cell, and will be automatically evaluated if you answer "Yes" to the initialization request.  If you do not answer "Yes" to this request, you must evaluate it before creating the graphs in the second section below.

Needs["Graphics`Animation`"]

Eliminating some unnecessary warning messages.  (Should be skipped)

Off[General::"spell"];
Off[General::"spell1"];

A first example (which you can easily modify)

The equation and the recurrence relation

Here we try to solve the differential equation

(2 + (x^2)) y" -  x y' + 4 y = 0.

0 is an ordinary point for this equation,  
    so we assume a solution of the form

y (x) = (∑_(n = 0)^∞) 
              (Subscript[a, n]) (x^n).

The differential equation gives us

(2 + (x^2)) (∑_(n = 2)^∞) 
              n (n - 1) (Subscript[a, n]) (x^(n - 2)) - 
        x (∑_(n = 1)^∞) n (Subscript[a, n]) (x^(n - 1)) + 
        4 (∑_(n = 0)^∞) (Subscript[a, n]) (x^n) = 0

which (after index shifts) reduces to

(∑_(n = 0)^∞)(2 (n + 2) (n + 1) 
                (Subscript[a, n + 2]) + n (n - 1) (Subscript[a, n]) - n (Subscript[a, n]) + 4 (Subscript[a, n])) 
          (x^n) = 0.

The recurrence relation is then

2 (n + 2) (n + 1) (Subscript[a, n + 2]) + n (n - 1) (Subscript[a, n]) - 
        n (Subscript[a, n]) + 4 (Subscript[a, n]) = 0.

or equivalently

(Subscript[a, n + 2]) = ((-((n^2) - 2  n + 
                  4))/(2 (n + 2) (n + 
                  1))) (Subscript[a, n]).

We can generate as many coefficients as we wish.  To get up to the coefficient of (x^k) (you can change k), use

Clear[a, k, n, x]
k=10;
a[0] = a0;
a[1] = a1;
Prepend[Prepend[Table[a[
                        n + 2] = ((-((n^2) - 2  n + 4) a[n])/(2 (n + 2) (n + 1))),{n, 0, k - 2}],a[1]],a[0]].Table[(x^n), {n, 0, k}]

I don't think it will be too easy to get the general term in this example.

An initial value problem with animated graphs

Let's assume we have some initial conditions, say y(0) = -1 and y'(0) = 2.  This tells us that a0 = -1 and a1=2.  Once again we can generate polynomial approximations of any degree (within computational reason).

Clear[a, k, n, x]
k=10;
a[0]=-1;
a[1]=2;
Prepend[Prepend[Table[a[
                        n + 2] = ((-((n^2) - 2  n + 4) a[n])/(2 (n + 2) (n + 1))),{n, 0, k - 2}],a[1]],a[0]].Table[(x^n), {n, 0, k}]

We can animate the graphs of polynomial approximations of increasing degree using the cell below.  To see the animation, evaluate the cell, select the entire collection of output graphs by clicking on the bracket containing all of them, and choose Animate Selected Graphics from the Cell menu item.  You can control the speed and direction of the animation using the VCR-like controls at the lower left corner of the notebook window.

m is the maximal degree polynomial to be plotted, a[0] and a[1] come from the initial conditions, xmin and xmax are the bounds on the x-axis, and ymin and ymax are the bounds on the y-axis.  A numerical solution is generated and graphed in black, while the polynomial approximation is graphed in red.

Clear[a, k, x, y, n, m, xmin, xmax, ymin, ymax, solution, 
      solPlot]
m=10;
a[0]=-1;
a[1]=2;
xmin=-2;
xmax=2;
ymin=-3;
ymax=3;
solution = NDSolve[{(2+(x^2)) y''[x] -  x y'[x] + 4 y[x]== 0,y[0]==a[0],y'[0]==a[1]},y,{x, xmin, xmax}];solPlot = 
        Plot[Evaluate[y[x]  /.  solution], {x, xmin, xmax}, 
          PlotRange -> {ymin, ymax}, DisplayFunction -> Identity];
Animate[Show[solPlot,Plot[Prepend[Prepend[Table[a[
                                      n + 2] = ((-((n^2) - 2  n + 4) a[n])/(2 (n + 2) (n + 1))),{n, 0, k - 2}],a[1]],a[0]].Table[(x^n), {n, 0, k}],{x, xmin, xmax},PlotRange -> {ymin, ymax},PlotStyle -> {RGBColor[1, 0, 0]},DisplayFunction -> Identity],PlotLabel -> "degree=" <> ToString[k],DisplayFunction -> $DisplayFunction],{k, 1, m, 1}];

Legendre's equation

The equation and the recurrence relation

Legendre's equation is the second order differential equation

(1 - (x^2)) y" -  2 x y' + p (p+1) y = 0.

0 is an ordinary point for this equation,  
    so we assume a solution of the form

y (x) = (∑_(n = 0)^∞) 
              (Subscript[a, n]) (x^n).

The differential equation gives us

(1 - (x^2)) (∑_(n = 2)^∞) 
              n (n - 1) (Subscript[a, n]) (x^(n - 2)) - 
        2  x (∑_(n = 1)^∞) n (Subscript[a, n]) (x^(n - 1)) + 
        p (p + 1) (∑_(n = 0)^∞) (Subscript[a, n]) (x^n) = 
      0

which (after index shifts) reduces to

(∑_(n = 0)^∞)((n + 2) (n + 1) 
                (Subscript[a, n + 2]) - n (n - 1) (Subscript[a, n]) - 2  n (Subscript[a, n]) + 
              p (p + 1) (Subscript[a, n])) (x^n) = 0.

The recurrence relation is then

(n + 2) (n + 1) (Subscript[a, n + 2]) - n (n - 1) (Subscript[a, n]) - 
        2  n (Subscript[a, n]) + p (p + 1) (Subscript[a, n]) = 0.

or equivalently

(Subscript[a, n + 2]) = ((-(p - n) (p + n + 
                  1))/((n + 2) (n + 
                  1))) (Subscript[a, n]).

We can generate as many coefficients as we wish.  To get up to the coefficient of (x^k) (you can change k), use

Clear[p, a, k, n, x]
k=10;
a[0] = a0;
a[1] = a1;
Prepend[
        Prepend[Table[
            a[n + 2] = ((-(p - n) (p + n + 1) a[n])/((n + 2) (n + 1))), {n, 0, k - 2}], a[1]], a[0]] . 
      Table[(x^n), {n, 0, k}]

Polynomial solutions and series solutions with animated graphs

We can get two linearly independent solutions by setting y(0)=a[0]=1 and y'(0)=a[1]=0 for one solution and y(0)=a[0]=0 and y'(0)=a[1]=1 for another one.  Notice that one of the two solutions we obtain this way is a terminating series (and therefore a polynomial) if p is a non-negative integer.  For even p it is the first of these that terminates, while for odd p it is the second.  Such solutions are called the Legendre polynomials (once they have been scaled by multiplying by a constant for reasons you may, or may not, see later).  You can obtain them by varying p, a[0], and a[1] in the cell below.  The default values produce the second Legendre polynomial if we set a[0] = -(1/2) instead of 1.  (Notice that although we set k = 10, only terms up to (x^2) appear because the series terminates.)

Clear[p, a, k, n, x, xmin, xmax, solution]
p=2;
k=10;
a[0]=-(1/2);
a[1]=0;
xmin=-1;
xmax=1;
solution = 
      Prepend[Prepend[
            Table[a[n + 
                    2] = ((-(p - n) (p + n + 1) a[n])/((n + 2) (n + 1))), {n, 0, k - 2}], a[1]], a[0]] . 
        Table[(x^n), {n, 0, k}]
Plot[
        solution, {x, xmin, xmax}];

Mathematica, of course, has the Legendre polynomials as built-in functions, so you can compare with the cell above and see what the value of a[0] needs to be to get the actual Legendre polynomial .

LegendreP[p, x]

We can animate the graphs of polynomial approximations of increasing degree using the cell below.  This is more interesting if the solution is one corresponding to an infinite series rather than a polynomial (whose graph we can produce in the cell above anyway).  To see the animation, evaluate the cell, select the entire collection of output graphs by clicking on the bracket containing all of them, and choose Animate Selected Graphics from the Cell menu item.  You can control the speed and direction of the animation using the VCR-like controls at the lower left corner of the notebook window.

m is the maximal degree polynomial to be plotted, a[0] and a[1] come from the initial conditions, xmin and xmax are the bounds on the x-axis, and ymin and ymax are the bounds on the y-axis.  A numerical solution is generated and graphed in black, while the polynomial approximation is graphed in red.  The default values produce the "Legendre function of the second kind of order p = 2" (i.e. the infinite series solution) if we set a[1] = -2 instead of 1.

Clear[p, a, k, x, y, n, m, xmin, xmax, ymin, ymax, solution, 
      solPlot]
p=2;
m=10;
a[0]=0;
a[1]=-2;
xmin=-1;
xmax=1;
ymin=-1;
ymax=1;
solution = NDSolve[{(1-(x^2)) y''[x] -  2  x y'[x] + p (p + 1) y[x]== 0,y[0] == a[0],y'[0] == a[1]},y,{x, xmin, xmax}];solPlot = 
        Plot[Evaluate[y[x]  /.  solution], {x, xmin, xmax}, 
          PlotRange -> {ymin, ymax}, DisplayFunction -> Identity];
\(Animate[
        Show[solPlot, 
          Plot[Prepend[
                Prepend[
                  Table[a[
                        n + 2] = \(\(-\((p - n)\)\) \((p + n + 1)\) a[n]\)\/\(\((n + 2)\) \((n + 1)\)\), {n, 0, k - 2}], a[1]], a[0]] . 
              Table[x\^n, {n, 0, k}], {x, xmin, xmax}, 
            PlotRange -> {ymin, ymax}, 
            PlotStyle -> {RGBColor[1, 0, 0]}, 
            DisplayFunction -> Identity], 
          PlotLabel -> "degree=" <> ToString[k], 
          DisplayFunction -> $DisplayFunction], {k, 1, m, 
          1}];\)

Once again, Mathematica has solutions of this type as built-in functions.  We chose a[1]=-2 in the cell above to make them match for p = 2.  The next cell compares the solution graph we got to the graph of LegendreQ[p,x].  The first graph is that of the function LegendreQ[p,x] and the second overlays the solution plot from the previous cell on top of it.  Make sure that cell was evaluated before evaluating this one.

Clear[x, legendrePlot]
legendrePlot = 
        Plot[LegendreQ[p, x], {x, xmin, xmax}, 
          PlotStyle -> RGBColor[0, 1, 0], 
          PlotRange -> {ymin, ymax}];
Show[solPlot, legendrePlot, 
        DisplayFunction -> $DisplayFunction];

Bessel's equation

Bessel's equation is the second order differential equation

(x^2) y" + x y' + ((x^2) - (p^2)) y = 0.

\(TextForm\`We\ will\ assume\ p\ is\ positive . \ \ 0\ is\ a\ regular\ singular\ point\ for\ this\ equation, \ 
    so\ we\ assume\ a\ solution\ of\ the\ form\)

y (x) = (∑_(n = 0)^∞) 
              (Subscript[a, n]) (x^(n + r)).

The differential equation gives us

(x^2) (∑_(n = 0)^∞)(n + r) (n + 
                    r - 1) (Subscript[a, n]) (x^(n + r - 2))+x (∑_(n = 1)^∞)(n + r) (Subscript[a, n]) 
                (x^(n + r - 1))+((x^2)-(p^2))(∑_(n = 0)^∞) (Subscript[a, n]) 
                (x^(n + r))=0

which (after index shifts) becomes

(r (r - 1)+r-(p^2))(Subscript[a, 0])+((r + 1) r+r+1-(p^2))(Subscript[a, 1])+(∑_(n = 2)^∞)((n + r) (n + r - 1) (Subscript[a, n])+(n + r) (Subscript[a, n])+ (Subscript[a, n - 2])-(p^2)(Subscript[a, n]))(x^(n + r))=0.

Simplifying a bit we get

((r^2)-(p^2))(Subscript[a, 0]) (x^r)+((((r + 1))^2)-(p^2))(Subscript[a, 1])(x^(1 + r))+(∑_(n = 2)^∞)(((((n + r))^2)-(p^2))(Subscript[a, n])+ (Subscript[a, n - 2]))(x^(n + r))=0.

The indicial equation is

(r^2)-(p^2)=0.

so r = p or r = -p.

We also know that the term ((((r+1))^2)-(p^2)) (Subscript[a, 1]) has to be zero.  If p ≠ (1/2) or r = p, then (Subscript[a, 1]) has to be zero.  If p = (1/2) and r = -(1/2), then we have a free choice of (Subscript[a, 1]), which we can set equal to zero.

The recurrence relation is

((((n + r))^2)-(p^2))(Subscript[a, n])+ (Subscript[a, n - 2])=0.

or equivalently

(Subscript[a, n])=((-1)/((((n + r))^2)-(p^2)))(Subscript[a, n - 2]).

as long as (((n+r))^2) - (p^2) ≠ 0.  The only time this can be zero is if r = -p and 0 = (((n-p))^2) - (p^2) = (n^2) - 2 n p = n (n - 2p).  Since this is for n ≥ 2, we are fine unless 2p is an integer.

In fact, we don't even have a problem if 2p is an odd integer.  In this case with r = -p we only run into trouble finding (Subscript[a, 2p]).  Since all of the odd coefficients up to that point will have to be 0 (since (Subscript[a, 1])= 0), we can simply set this one to be 0 as well, and the recurrence relation is still satisfied.  As a result, all the odd coefficients will be 0, and the only problems will be when p is an integer.

Corresponding to the larger root (r = p) we always get a solution to Bessel's equation.  In this case (((n+r))^2) - (p^2) simplifies to n(n + 2p), so the recurrence relation becomes

(Subscript[a, n]) = ((-1)/(n (n + 
                  2  p))) (Subscript[a, n - 2]).

We can generate a few terms using Mathematica.

Clear[p, a, a0, k, n, x]
k=10;
a[0] = a0;
a[1] = 0;
Prepend[
        Prepend[Table[a[n] = ((-a[n - 2])/(n (n + 2  p))), {n, 2, k}], 
          a[1]], a[0]] . Table[(x^(n + p)), {n, 0, k}]

What may not be so clear from this is a general formula for the integers that are at the beginning of each denominator.  From the recurrence relation it is not too hard to see that in the coefficient of (x^(2 n + p)) that integer is the product of the first n even integers.  (This comes from the leading n in the denominator of our recurrence relation.)  This, in turn, can be expressed as (2^n)n! by factoring out the 2s in the product of even integers.  We can verify this using Mathematica.

Clear[n]
Table[(2^n) n!, {n, 0, ⌊(k/2)⌋}]

So a general term of the solution has the form

(((((Subscript[a, 0])(-1))^n)(x^(2n+p)))/((2^n)n!(2+2p)(4+2p)...(2n+2p))).

We can factor 2 out of each factor of the denominator that involves p (n such factors, producing a (2^n) we can combine with the one already in the denominator) to get

(((((Subscript[a, 0])(-1))^n)(x^(2n+p)))/((2^(2n))n!(1+p)(2+p)...(n+p))).

Remember that for positive integers n and p we know that

p!(1+p)(2+p)...(n+p) = (n+p)!.

We can use the same notation even if p is not an integer, in which case

(1+p)(2+p)...(n+p) = (((n+p)!)/(p!)).

Mathematica actually knows this notation even for non-integers p, as you can see by the following example.

Clear[n, p, k]
p=(1/3);
n=5;
N[(∏_(k = 1)^n)(k + 
            p)]
N[(((n + p)!)/(p!))]

(The definition of p! for non-integers p involves a function called the gamma function, denoted by Γ(x) and defined by

Γ(x) = (∫_0^∞)(t^(x-1))(e^(-t))dt

for x > 0.  This definition can be extended to all x except non-positive integers.  It can be shown that Γ(p+1) = p! for positive integers, so we can use this to define p! for all p which are not negative integers.)

After all this, we can rewrite the general term as

(((((Subscript[a, 0])(-1))^n)(x^(2n+p))p!)/((2^(2n))n!(n+p)!)).

We can freely choose (Subscript[a, 0])≠0.  If we make (Subscript[a, 0]) = (1/((2^p)p!)), we can rewrite this as

(((((-1))^n)(x^(2n+p)))/((2^(2n+p))n!(n+p)!))

and finally as

(((((-1))^n)((((x/2)))^(2n+p)))/(n!(n+p)!)).

This means we have as one solution of Bessel's equation the function

(Subscript[J, p])(x) = (∑_(n=0)^∞)(((((-1))^n)((((x/2)))^(2n+p)))/(n!(n+p)!)).

This function is called the Bessel function of the first kind of order p.

We are still faced with finding a second linearly independent solution.  We know from our solution method (method of Frobenius) that all the above work will go through toward finding a second solution as long as the roots of the indicial equation (which were p and -p) do not differ by an integer.  It turns out that there are two different cases to consider.

p is not an integer

We can animate the convergence of the series to (Subscript[J, p])(x) using the cells below.  For each graph we keep terms up to ((((x/2)))^(2  k+p)).  Notice that Mathematica refers to the Bessel function of the first kind of order p by BesselJ[p,x].  The graphs of the approximating partial sums are in red and the graph of (Subscript[J, p])(x) is in blue.

Clear[p, k, x, y, n, m, xmin, xmax, ymin, ymax, 
      solPlot1]
p=(1/2);
m=10;
xmin=0;
xmax=10;
ymin=-2;
ymax=2;solPlot1 = 
        Plot[BesselJ[p, x], {x, xmin, xmax}, PlotRange -> {ymin, ymax}, 
          PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
\(Animate[
        Show[solPlot1, 
          Plot[∑\+\(n = 0\)\%k\(\(\((\(-1\))\)\^n\) \((x\/2)\)\^\((2  n + p)\)\)\/\(\(n!\) \(\((n + p)\)!\)\), {x, xmin, xmax}, 
            PlotRange -> {ymin, ymax}, 
            PlotStyle -> {RGBColor[1, 0, 0]}, 
            DisplayFunction -> Identity], 
          PlotLabel -> 
            ToString[2  k + p] <> " is the highest power", 
          DisplayFunction -> $DisplayFunction], {k, 1, m, 
          1}];\)

In this situation everything we did in the section above goes through for the smaller root (r = -p), and so (Subscript[J, -p])(x) is a second linearly independent solution.  Notice that the initial term of  (Subscript[J, -p])(x) is (1/((-p)!))((((x/2)))^(-p)), which is unbounded at 0, so (Subscript[J, p])(x) and (Subscript[J, -p])(x) are clearly independent.  Here's the animation for the convergence to the second solution.  Make sure the cell above has been evaluated during the current Mathematica session before evaluating this one.  The graphs of the approximating partial sums are in red and the graph of (Subscript[J, -p])(x) is in green.

Clear[r, k, x, y, n, m, 
      solPlot2]
r = -p;
m=10;
solPlot2 = 
        Plot[BesselJ[r, x], {x, xmin, xmax}, PlotRange -> {ymin, ymax}, 
          PlotStyle -> RGBColor[0, 1, 0], 
          DisplayFunction -> Identity];
\(Animate[
        Show[solPlot2, 
          Plot[∑\+\(n = 0\)\%k\(\(\((\(-1\))\)\^n\) \((x\/2)\)\^\((2  n + r)\)\)\/\(\(n!\) \(\((n + r)\)!\)\), {x, xmin, xmax}, 
            PlotRange -> {ymin, ymax}, 
            PlotStyle -> {RGBColor[1, 0, 0]}, 
            DisplayFunction -> Identity], 
          PlotLabel -> 
            ToString[2  k + r] <> " is the highest power", 
          DisplayFunction -> $DisplayFunction], {k, 1, m, 
          1}];\)

Provided the two input cells above have been evaluated, we can plot the two solutions together on the same set of axes.

Show[solPlot1, solPlot2, 
        DisplayFunction -> $DisplayFunction];

By the way, a (possibly different) second linearly independent solution for Bessel's equation of order p (instead of (Subscript[J, -p])(x)) is the Bessel function of the second kind of order p denoted by (Subscript[Y, p])(x).  Mathematica uses BesselY[p,x] for this function.  When p is not an integer, (Subscript[Y, p])(x) = (((Subscript[J, p])(x) cos(p π)- (Subscript[J, -p])(x))/(sin(p π))) (strange, but there are reasons for this).  Notice that if p is an odd multiple of (1/2) then (Subscript[Y, p])(x) = ± (Subscript[J, -p])(x).  We can see the equality using the cell below.  The first graph is (Subscript[Y, p])(x), the second is (((Subscript[J, p])(x) cos(p π)- (Subscript[J, -p])(x))/(sin(p π))), and the third is the overlay of the two.

Clear[p]
p=(1/2);
Show[Plot[BesselY[p, x], {x, xmin, xmax}], 
        Plot[((BesselJ[p, x] Cos[p π] - BesselJ[-p, x])/Sin)[p π], {x, xmin, xmax}]];

One last interesting note.  When p = (1/2), our general term

(((((Subscript[a, 0])(-1))^n)(x^(2n+p)))/((2^n)n!(2+2p)(4+2p)...(2n+2p))).

with which we began our manipulations is

(((((Subscript[a, 0])(-1))^n)(x^(2n+(1/2))))/((2^n)n!(3)(5)...(2n+1))),

which can be rewritten as

(Subscript[a, 0])(x^(-(1/2)))(((((-1))^n)(x^(2n+1)))/((2n+1)!)),

and therefore the series is that of

(Subscript[a, 0])(x^(-(1/2))) sin(x).

If we choose (Subscript[a, 0]) = (Sqrt[(2/π)]) we get (Subscript[J, (1/2)])(x).

When p = -(1/2), our general term is

(((((Subscript[a, 0])(-1))^n)(x^(2n-(1/2))))/((2^n)n!(1)(3)...(2n-1))),

which can be rewritten as

(Subscript[a, 0])(x^(-(1/2)))(((((-1))^n)(x^(2n)))/((2n)!)),

and therefore the series is that of

(Subscript[a, 0])(x^(-(1/2))) cos(x).

Again if we choose (Subscript[a, 0]) = (Sqrt[(2/π)]) we get (Subscript[J, -(1/2)])(x).

(In fact, for large values of x any Bessel function (not just order (1/2)) has limiting behavior like some linear combination of (x^(-(1/2))) sin(x) and (x^(-(1/2))) cos(x).  We will not investigate this here.)

p is an integer

Once again corresponding to the larger root (r = p) we have one solution (Subscript[J, p])(x).  We can animate the convergence of the series to (Subscript[J, p])(x) using the cells below.  For each graph we keep terms up to ((((x/2)))^(2  k+p)).  The graphs of the approximating partial sums are in red and the graph of (Subscript[J, p])(x) is in blue.

Clear[p, k, x, y, n, m, xmin, xmax, ymin, ymax, 
      solPlot1]
p=0;
m=10;
xmin=0;
xmax=10;
ymin=-1;
ymax=1;solPlot1 = 
        Plot[BesselJ[p, x], {x, xmin, xmax}, PlotRange -> {ymin, ymax}, 
          PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
\(Animate[
        Show[solPlot1, 
          Plot[∑\+\(n = 0\)\%k\(\(\((\(-1\))\)\^n\) \((x\/2)\)\^\((2  n + p)\)\)\/\(\(n!\) \(\((n + p)\)!\)\), {x, xmin, xmax}, 
            PlotRange -> {ymin, ymax}, 
            PlotStyle -> {RGBColor[1, 0, 0]}, 
            DisplayFunction -> Identity], 
          PlotLabel -> 
            ToString[2  k + p] <> " is the highest power", 
          DisplayFunction -> $DisplayFunction], {k, 1, m, 
          1}];\)

Unfortunately, when p is an integer (Subscript[J, -p])(x) turns out to be (((-1))^p) (Subscript[J, p])(x), and so does not give us a second solution.  We can see this relationship using the cell below.  The first graph is that of (Subscript[J, -p])(x), the second is of (((-1))^p) (Subscript[J, p])(x), and the third is the overlay of the two.

Clear[p, xmin, xmax, ymin, ymax]
p=1;
xmin=0;
xmax=10;
ymin=-1;
ymax=1;
Show[Plot[BesselJ[-p, x], {x, xmin, xmax}], 
        Plot[(((-1))^p) BesselJ[p, x], {x, xmin, 
            xmax}]];

(Actually, from our work above it is not even clear that we can define (Subscript[J, -p])(x) by using the series we created, as the (n-p)! in the denominator is not defined for n < p.  It turns out that this does not really pose a problem as the gamma function approaches either -∞ or ∞ using one sided limits near non-positive integers, and so the reciprocal approaches 0 there.  It doesn't matter for our purposes - we still do not get a second linearly independent solution this way.)

A second linearly independent solution can be found in several different ways, but we will not do any of them here.  Mathematica does know a second linearly independent solution, a Bessel function of the second kind of order p, again called (Subscript[Y, p])(x).  The cell below draws the graphs of (Subscript[J, p])(x) (in blue) and (Subscript[Y, p])(x) (in green).

Clear[p, xmin, xmax, ymin, ymax, solPlot1, solPlot2]
p=1;
xmin=0;
xmax=10;
ymin=-1;
ymax=1;
solPlot1 = 
        Plot[BesselJ[p, x], {x, xmin, xmax}, PlotRange -> {ymin, ymax}, 
          PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
solPlot2 = 
        Plot[BesselY[p, x], {x, xmin, xmax}, PlotRange -> {ymin, ymax}, 
          PlotStyle -> RGBColor[0, 1, 0], 
          DisplayFunction -> Identity];
Show[solPlot1, 
        solPlot2, DisplayFunction -> $DisplayFunction];