# -*- coding: utf-8 -*-
"""
Created on Thu Oct 16 07:21:55 2014

@author: schildi
"""

#Imports 
from matplotlib.pyplot import show, plot, legend, xlabel, ylabel
from numpy import linspace, zeros, r_, tanh, array, sqrt
from scipy.optimize import fsolve


#Problem definition, y1 has to be adjusted to the value of the 

xi=1.0
F=lambda y: array([y[1], 1.0*(y[0]**3-y[0])/(xi**2)])
y0=0
y1=1.0/(sqrt(2)*xi)


#Parameters for programm

Integration_Depth = 4
Sample_Points=250
X=linspace(0,Integration_Depth,Sample_Points)
h=X[1]-X[0]
Y=zeros((Sample_Points,2))
Y[0,0]=y0
Y[0,1]=y1

#use explicit Euler (could be better but that in case if better needed)

for i in r_[0:Sample_Points-1]:
    Y[i+1,:]=Y[i]+h*F(Y[i,:])
    print Y[i+1,:]

    
#use implicit midpoint rule because of critical point at f=1
yIMP=zeros((Sample_Points,2))
yIMP[0,0]=y0
yIMP[0,1]=y1
for i in r_[0:Sample_Points-1]:
    ftilde = lambda y: yIMP[i]+h*F((y+yIMP[i,:])/2.)-y
    yIMP[i+1,:]=fsolve(ftilde, yIMP[i,:])


#Extract solution

Ysol=Y[:,0]
YsolIMP=yIMP[:,0]

#Plot solution and the right function
figure(1)
plot(X,Ysol, 'r', label='Numerical Solution EE')
plot(X,YsolIMP, 'b', label='Numerical Solution IMP')
xlabel('length / xi')
ylabel('f')
X1=X+0.5*h
plot(X1,tanh(X1/(sqrt(2)*xi)),'g', label='exact solution out of lecture')
legend(loc=0)
show()

#Plot the difference between right solution and the numerical result
figure(2)
YErr=Ysol-tanh(X/(sqrt(2)*xi))
YErrIMP=YsolIMP-tanh(X/(sqrt(2)*xi)) 
#plot(X,YErr, label='Error EE')
plot(X,YErrIMP, label='Error IMP')
xlabel('length / xi')
ylabel('numerical result - exact result')
legend(loc=0)
show()