Untitled

import scipy as S
import numpy as N
import numpy as N
from numpy import allclose, arange, eye, linalg, ones, zeros
from scipy.sparse.linalg import dsolve
from scipy import sparse, matrix
from decimal import *
getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999999, Emax=999999999,
        capitals=1, flags=[], traps=[Overflow, DivisionByZero,
        InvalidOperation])
re=5000
n=100
u=zeros((n+1,1), dtype=float)
k=zeros((n+1,1), dtype=float)
unew=zeros((n+1,1), dtype=float)
iter=3000.0
float (4)
dy=1.0/n
dt=re*dy*dy

u[n]=1
b=1.0+dt/dy/dy/re
a=-dt/2.0/dy/dy/re
t=0.0
for i in range(1, n):
    t=t+dt
    for j in range (2, n):
        k[j]=(1.0-dt/dy/dy/re)*u[j]+dt/2.0/dy/dy/re*(u[j+1]+u[j-1])
    k[2]=k[2]-a*u[n]
    k[n]=k[n]-a*u[n+1]

    #przeskakujemy do linii 95 bo tam sie definiuje macierze do rozwiazania

    aa=zeros((n-1,n-1))  #linia 95
    bb=zeros((n-1),1)   #linia 95
    x=zeros((n-1),1)    #linia 95
    kk=zeros((n,1)) #linia 97
    for j in range(1,n-1)   #linia 99
        bb[j]=kk[j+1]
    for i in range (1,n-1): #linia 104
        for j in range(1,n-1):
            aa[i,j]=0.0;    #liina 106


        for j in range(1, n-1):
            u[j+1]=unew[j]
#koniec, teraz trzeba linie 73 odpalic