I figured out the following Python code - in Spyder - for a particular solution to the Heat-Diffusion equation, last week, which follows for the case when k = 4. I took the initial-boundary-condition to be u(x) = 4*cos(x), where u_x(x) = -4*sin(x). Please note that I had to cut-and-paste the final output solution: u(t, x) = 4*(sin(x + pi/2) + 3*sin(3*x + pi/2))*(-sinh(4*t) + cosh(4*t)) into the plot commands at the end, but this does seem to indicate suitable ways and means of ANALYTICALLY solving SECOND-ORDER PDEs reasonably successfully in Spyder. I hope that you find this program both interesting and potentially useful. I have also written codes for the wave-equation and the Schrodinger equation, both of which work well. The following code, with a Word document attached showing the plot itself, is:

from __future__ import division

from sympy import *

import numpy as np

import sympy as sy

import matplotlib.pyplot as plt

from sympy import init_printing

init_printing()

from sympy import Function, Indexed, Tuple, sqrt, dsolve, solve, Eq, integrate,

Derivative, sin, cos, symbols

from sympy.abc import k

from sympy import solve, Poly, Eq, Function, exp

from sympy.abc import x, y, z, a, b, c, k

f = Function('f')

from sympy import Indexed, IndexedBase, Tuple, sqrt

from IPython.display import display

from sympy import *

from sympy.abc import *

from sympy.plotting import plot

init_printing()

# The code below solves the Heat-Diffusion PDE, uxx = ut, by first separating

# it out using eigenfunctions,and then solving the associated ODEs for X(t) and

# T(t) to then yield the final solution for u(t, x), where u(t, x) = X(x)T(t)

import sympy

sympy.init_printing()

# First, we define our symbols and function:

from sympy import Function, pprint

from sympy import Function, Eq, pde_separate_mul, Derivative as D

from sympy.abc import x, t

u, X, T = map(Function, 'uXT')

# Separating out the PDE, we have

eq = Eq(D(u(x, t), x, 2), D(u(x, t), t, 1))

display(pde_separate_mul(eq, u(x, t), [X(x), T(t)]))

# The solution for X(x)

m, k, u_0, r_0, x = sympy.symbols('m, k, u_0, r_0, x')

X = sympy.Function('X')

sol = sympy.dsolve(sympy.Derivative(X(x), x, 2) + k * X(x))

sol

display(simplify(sol))

# And define our initial conditions and solve for them:

ics = [sympy.Eq(sol.args[1].subs(x, 0), 4*cos(x)),

sympy.Eq(sol.args[1].diff(x).subs(x, 0), -4*sin(x))]

ics

# Particular solution for general values of k

solved_ics = sympy.solve(ics)

solved_ics

display(simplify(solved_ics))

full_sol = sol.subs(solved_ics[0])

display(simplify(full_sol))

# Second part of PDE in terms of time, t

t, C1, C2, C3 = symbols("t C1 C2 C3")

x, y, z = symbols("x y z", cls = Function, Function = True)

import sympy as sm

import numpy as np

import matplotlib.pyplot as plt

sm.init_printing()

from sympy.abc import t, k

T = sm.Function('T')(t)

eq1 = sm.Eq(T.diff(t), -k*T)

sol3 = sm.dsolve(eq1)

sol3

# Initial conditions

t0 = sol3.args[1].subs({'t': 0})

T0 = sm.symbols('4')

eq_init = sm.Eq(T0, t0)

# Particular solution for general values of k

C1 = eq_init.args[1]

sm.solve(eq_init, C1)

init_solve = sm.solve(eq_init, C1)

full_sol4 = sol3.subs(C1, init_solve[0])

display(simplify(full_sol4))

display((simplify((full_sol).rhs*(full_sol4).rhs)).rewrite(sympy.sin))

# Particular solution for the original PDE for specific parameter values of k

from sympy import Symbol

k = Symbol('k')

f = (((full_sol).rhs*(full_sol4).rhs).rewrite(sympy.sin))

display((simplify(f.subs({k:4}))).rewrite(sympy.sin))

print((simplify(f.subs({k:4}))).rewrite(sympy.sin))

# Plot 1

from sympy import Function, pprint

import numpy as np

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

import random

from numpy import sqrt, pi, exp, sin, cos, sinh, cosh, arange

from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show

def fun(x, t):

# The next line requires an indent

return 4*(sin(x + pi/2) + 3*sin(3*x + pi/2))*(-sinh(4*t) + cosh(4*t))

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

x = t = np.arange(-2.0, 2.0, 0.05)

X, T = np.meshgrid(x, t)

zs = np.array([fun(x,t) for x,t in zip(np.ravel(X), np.ravel(T))])

Z = zs.reshape(X.shape)

ax.plot_surface(X, T, Z)

ax.set_xlabel('Distance, x')

ax.set_ylabel('Time, t')

ax.set_zlabel('u(t, x)')

plt.show()

For those who are interested, I have a similar Python solution here to the Schrodinger equation for the initial-boundary condition u(x) = 3cos(x), with derivative u_x(x) = -3sin(x) (although, you could try u(x) = 3 and u_x(x) = -3, depending on what is required). You should verify if the solution to this code is correct, so do let me know if there is anything amiss. A Word document containing the two plots that it generates is attached.

import numpy as np

import sympy as sy

import matplotlib.pyplot as plt

from sympy import init_printing

import cmath

init_printing()

from sympy import Function, Indexed, Tuple, sqrt, dsolve, solve, Eq, integrate, Derivative, sin, cos, symbols

from sympy.abc import k

from sympy import solve, Poly, Eq, Function, exp

from sympy.abc import x, y, z, a, b, c, f, E, V

f = Function('f')

from sympy import Indexed, IndexedBase, Tuple, sqrt

from IPython.display import display

from sympy import *

from sympy.abc import *

from sympy.plotting import plot

from sympy import trigsimp, cancel

init_printing()

# The code below solves the wave-equation PDE, uxx = utt, by first separating

# it out using eigenfunctions,and then solving the associated ODEs for X(t) and

# T(t) to then yield the final solution for u(t, x), where u(t, x) = X(x)T(t)

import sympy

sympy.init_printing()

# First, we define our symbols and function:

from sympy import Function, pprint

from sympy import Function, Eq, pde_separate_mul, Derivative as D

from sympy.abc import x, t, f

u, X, T = map(Function, 'uXT')

# Run code block in SymPy Live

eq = Eq(D(u(x, t), x, 2) + V*u(x, t), D(u(x, t), t, 2))

display(pde_separate_mul(eq, u(x, t), [X(x), T(t)]))

# The solution for X(x):

m, E, V, u_0, r_0, x = sympy.symbols('m, E, V, u_0, r_0, x')

X = sympy.Function('X')

sol = sympy.dsolve(sympy.Derivative(X(x), x, 2) + (V - E) * X(x))

sol

display(simplify(sol))

# The initial-boundary-conditions are:

ics = [sympy.Eq(sol.args[1].subs(x, 0), 3*cos(x)),

sympy.Eq(sol.args[1].diff(x).subs(x, 0), -3*sin(x))]

ics

# Particular solution for X(x) for general values of k

solved_ics = sympy.solve(ics)

solved_ics

display(simplify(solved_ics))

full_sol = sol.subs(solved_ics[0])

display(simplify(full_sol))

# Second part of PDE in terms of time, t

t, C1, C2, C3 = symbols("t C1 C2 C3")

x, y, z = symbols("x y z", cls = Function, Function = True)

import sympy as sm

import numpy as np

import matplotlib.pyplot as plt

sm.init_printing()

from sympy.abc import t, k

T = sm.Function('T')(t)

eq1 = sm.Eq(T.diff(t), - E * T)

sol2 = sm.dsolve(eq1)

sol2

display(sol2)

# Initial conditions

t0 = sol2.args[1].subs({'t': 0})

T0 = sm.symbols('3')

eq_init = sm.Eq(T0, t0)

# Particular solution for general values of k

C1 = eq_init.args[1]

sm.solve(eq_init, C1)

init_solve = sm.solve(eq_init, C1)

full_sol3 = sol2.subs(C1, init_solve[0])

display((simplify(full_sol3)).rewrite(sympy.sin))

display((simplify((full_sol).rhs*(full_sol3).rhs)).rewrite(sympy.sin))

# Particular solution for the original PDE for specific parameter values of k

# Full particular solution for given values of E and V

from sympy import Symbol

f = simplify((((full_sol).rhs*(full_sol3).rhs)).rewrite(sympy.sin))

f

display((simplify(f.subs({E: 4, V: 1}))).rewrite(sympy.sin))

print((simplify(f.subs({E: 4, V: 1}))).rewrite(sympy.sin))

from sympy import Symbol

f = simplify((((full_sol).rhs*(full_sol3).rhs)).rewrite(sympy.sin))

f

display((simplify(f.subs({E: 1, V: 4}))).rewrite(sympy.sin))

print((simplify(f.subs({E: 1, V: 4}))).rewrite(sympy.sin))

# Plot of each respective function:

# Plot 1

from sympy import Function, pprint

import numpy as np

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

import random

from numpy import sqrt, pi, exp, sin, cos, sinh, cosh, arange

from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show

def fun(x, t):

# The next line needs to be indented

return -3*(sqrt(3)*sin(x)*sinh(sqrt(3)*x) - 3*sin(x + pi/2)*cosh(sqrt(3)*x))*(-sinh(4*t) + cosh(4*t))

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

x = t = np.arange(-2.0, 2.0, 0.05)

X, T = np.meshgrid(x, t)

zs = np.array([fun(x,t) for x,t in zip(np.ravel(X), np.ravel(T))])

Z = zs.reshape(X.shape)

ax.plot_surface(X, T, Z)

# Plot 2

ax.set_xlabel('Distance, x')

ax.set_ylabel('Time, t')

ax.set_zlabel('u(t, x)')

plt.show()

from sympy import Function, pprint

import numpy as np

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

import random

from numpy import sqrt, pi, exp, sin, cos, sinh, cosh, arange

from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show

def fun(x, t):

# The next line needs to be indented

return 3*(-sqrt(3)*sin(x)*sin(sqrt(3)*x) + 3*sin(x + pi/2)*sin(sqrt(3)*x + pi/2))*(-sinh(t) + cosh(t))

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

x = t = np.arange(-2.0, 2.0, 0.05)

X, T = np.meshgrid(x, t)

zs = np.array([fun(x,t) for x,t in zip(np.ravel(X), np.ravel(T))])

Z = zs.reshape(X.shape)

ax.plot_surface(X, T, Z)

ax.set_xlabel('Distance, x')

ax.set_ylabel('Time, t')

ax.set_zlabel('u(t, x)')

plt.show()

And...here is a Python solution for the wave-equation with arbitrary initial-boundary-conditions u(x) = 3 and derivative u_x(x) = -3 and u(t) = 2 and derivative u_t(t) = -2 Again, you will find a plot in the Word document attached.

import numpy as np

import sympy as sy

import matplotlib.pyplot as plt

from sympy import init_printing

init_printing()

from sympy import Function, Indexed, Tuple, sqrt, dsolve, solve, Eq, integrate,

Derivative, sin, cos, symbols

from sympy.abc import k

from sympy import solve, Poly, Eq, Function, exp

from sympy.abc import x, y, z, a, b, c, k

f = Function('f')

from sympy import Indexed, IndexedBase, Tuple, sqrt

from IPython.display import display

from sympy import *

from sympy.abc import *

from sympy.plotting import plot

init_printing()

# The code below solves the wave-equation PDE, uxx = utt, by first separating

# it out using eigenfunctions,and then solving the associated ODEs for X(t) and

# T(t) to then yield the final solution for u(t, x), where u(t, x) = X(x)T(t)

import sympy

sympy.init_printing()

# First, we define our symbols and function:

from sympy import Function, pprint

from sympy import Function, Eq, pde_separate_mul, Derivative as D

from sympy.abc import x, t

u, X, T = map(Function, 'uXT')

# Run code block in SymPy Live

eq = Eq(D(u(x, t), x, 2), D(u(x, t), t, 2))

display(pde_separate_mul(eq, u(x, t), [X(x), T(t)]))

# The solution for X(x):

m, k, u_0, r_0, x = sympy.symbols('m, k, u_0, r_0, x')

X = sympy.Function('X')

sol = sympy.dsolve(sympy.Derivative(X(x), x, 2) + k * X(x))

sol

display(sol)

# The initial conditions are:

ics = [sympy.Eq(sol.args[1].subs(x, 0), 3), sympy.Eq(sol.args[1].diff(x).subs(x, 0), -3)]

ics

# Particular solution for X(x) for general values of k

solved_ics = sympy.solve(ics)

solved_ics

display(solved_ics)

full_sol = sol.subs(solved_ics[0])

display(full_sol)

# Second part of PDE in terms of time, t

m, k, v_0, r_0, t = sympy.symbols('m, k, v_0, r_0, t')

T = sympy.Function('T')

sol2 = sympy.dsolve(sympy.Derivative(T(t), t, 1) + k * T(t))

sol2

display(sol2)

# The initial-boundary-conditions are:

ics = [sympy.Eq(sol2.args[1].subs(t, 0), 2), sympy.Eq(sol2.args[1].diff(t).subs(t, 0), -2)]

ics

# Particular solution for T(t) general values of k

solved_ics = sympy.solve(ics)

solved_ics

display(solved_ics)

full_sol2 = sol2.subs(solved_ics[0])

full_sol2

display(full_sol2)

# Full particular solution for general values of k

display(((full_sol).rhs*(full_sol2).rhs).rewrite(sympy.sin))

print(((full_sol).rhs*(full_sol2).rhs).rewrite(sympy.sin))

# Full particular solution for a given value of k

from sympy import Symbol

k = Symbol('k')

f = (((full_sol).rhs*(full_sol2).rhs).rewrite(sympy.sin))

display((simplify(f.subs({k:100}))).rewrite(sympy.sin))

print((simplify(f.subs({k:100}))).rewrite(sympy.sin))

# Plot 1

from sympy import Function, pprint

import numpy as np

from mpl_toolkits.mplot3d import Axes3D

import matplotlib.pyplot as plt

import random

from numpy import sqrt, pi, exp, sin, cos, sinh, cosh, arange

from pylab import meshgrid,cm,imshow,contour,clabel,colorbar,axis,title,show

def fun(x, t):

# The next line requires an indent in front of the word 'return':

return 3*(sin(10*x) - 10*sin(10*x + pi/2))*(sinh(t) - cosh(t))/5

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

x = t = np.arange(-0.75, 0.75, 0.05)

X, T = np.meshgrid(x, t)

zs = np.array([fun(x,t) for x,t in zip(np.ravel(X), np.ravel(T))])

Z = zs.reshape(X.shape)

ax.plot_surface(X, T, Z)

ax.set_xlabel('Distance, x')

ax.set_ylabel('Time, t')

ax.set_zlabel('u(t, x)')

plt.show()