# Solving second-order ODE by the finite difference method in Python

Following code solves this second order linear ordinary differential equation $$ y''+7y=8\cos(4x)+\sin^{2}(2x), y(0)=\alpha, y(\pi/2)=\beta $$

by the finite differences method using just default libraries in `Python 3`

(tested with `Python 3.4`

). Linear system is solved by matrix factorization.

This snippet was used for NUM2 subject in FJFI, 2015 as a final project. Big thanks to my friend Vojta, who also participate.

**This is just for educational purposes and cannot be used for cheating**.

Here is the code (download):

```
from math import cos, sin, pi
# Settings
g1 = float(input("Value in y(0) = ")) # 3 ... y(0)
g2 = float(input("Value in y(pi/2) = ")) # 1.5 ... y(pi/2)
m = int(input("Number of steps: ")) # e.g. 1000
h = (pi/2)/m # krok site
net = [i*h for i in range(m+1)] # body site
# Initial
f = lambda x: 8 * cos(4*x) + (sin(2*x) ** 2) # f(x)
p = lambda x: -1 # p(x)
q = lambda x: 7 # q(x)
P = [p(i*h) for i in range(m+1)]
Q = [q(i*h) for i in range(m+1)]
F = [-f(i*h) for i in range(1,m)]
## Solving tridiagonal matrix by matrix factorization
A = [-(P[i+1]/h**2) for i in range(m-1)] # upper diagonal
C = [-(P[i+1]+P[i+2])/h**2 - Q[i+1] for i in range(m-1)] # diagonal
B = [-(P[i+2]/h**2) for i in range(m-1)] # lower diagonal
alphas = [0]
betas = [g1]
for i in range(m-1):
alphas.append(B[i] / (C[i] - alphas[i]*A[i]))
betas.append((betas[i]*A[i] + F[i]) / (C[i]-alphas[i]*A[i]))
u = [g2]
for i in range(m):
u.insert(0, alphas[m-i-1]*u[0]+betas[m-i-1])
## Write to file in gnuplot format
## In gnuplot just type: plot "out.txt"
with open("out.txt", "w", encoding="utf8") as outf:
outf.write("# x \t y\n")
for xova, yova in zip(net, u):
outs = "{0} {1}\n".format(xova, yova)
outf.write(outs)
```