import math as m
import matplotlib.pyplot as plt

u0 = 0.1
C = 30

def f1(s):
    return s**2

def f2(s):
    return (1-s)**2

def fi(s, c):
    a = mu(c) * ((s - 0.2) / 0.81)**3
    b = a + ((0.74 - s) / 0.715)**3
    return  a / b

def mu(c):
    return 20 * (1 + 5*c)

def a(c):
    return 0.1*c

tau = 0.0002
h = 0.01

N = 100
T = 1000
m = 0.2
U = 1
D = 0.0001

c = [[None for i in range(N)] for j in range(T)]
s = [[None for i in range(N)] for j in range(T)]

for i in range (N):
    s[0][i] = 0.74
    c[0][i] = 1

for j in range (T-1):
    s[j][0] = 0.2
    c[j][0] = 0
    c[j][N-1] = 1
    for i in range(1, N-1):
        s[j+1][i] = s[j][i] - (fi(s[j][i], c[j][i]) - fi(s[j][i-1], c[j][i-1]))/h*tau/m
        c1 = s[j][i]*c[j][i] + 0.1*c[j][i]/m + (-U*(fi(s[j][i], c[j][i]) - fi(s[j][i-1], c[j][i])*c[j][i-1]) + D*(c[j][i-1] - 2*c[j][i] + c[j][i+1])/h**2)*tau/m
        c[j+1][i] = c1/(s[j+1][i] + 0.1/m)

print(s[300])
print(c[300])