# -*- coding: utf-8 -*-
"""
Created on Wed Jan 21 17:19:19 201

Listing 12.1 : résolution de y''+ y = x, y(0) = y(1) = 0
@author: grivet
"""

from pylab import *

xmin = 0.0; xmax = 1.0
# ces valeurs ne sont pas des entiers
npt = int(input("nombre total de points: "))
nint = npt-1; nx = npt-2
M = zeros((nx,nx))
h = (xmax-xmin)/nint
sousdiag = ones(nx-1)
M = diag(sousdiag,-1)+diag(sousdiag,1)+(h*h-2)*diag(ones(nx),0)
x = linspace(xmin+h,xmax-h,nx)
c = h*h*x
y = solve(M,c)
yy = hstack([0,y,0]); xx = hstack([xmin,x,xmax])
plot(xx,yy,'ok')
zz = xx -(1/sin(1))*sin(xx)
plot(xx,zz,'-r')
title("Résolution d'un problème aux limites. \n Solution exacte(en rouge) et solution calculée (points noirs")