Hy, guys.
I write a subroutine which has to solve differential equation numerically. Here it is
subroutine Tridiagonal_Solv(n, a, b, alpha,beta)
implicit none
integer :: n, i
double precision, intent(in) :: a, b, alpha, beta
real :: h
double precision, dimension(0:5) :: x, r, u, v, w, z, s, t, y
h = Step(n, a, b)
x(1) = a
r(1) = alpha
do i = 2, n-1
x(i) = a+i*h
r(i) = 2*sin(x(i))
end do
r = beta
x = b
do i = 2, n-1
u(i) = -1/(h*h)
end do
u = 2/h
v(1) = -3/(2*h)
do i = 2, n-1
v(i) = 2/(h*h) + 1
end do
v = 3/(2*h)
w(1) = 2/h
do i = 2, n-1
w(i) = -1/(h*h)
end do
s(0) = 0
t(0) = 0
do i = 1, n
s(i) = -w(i)/(v(i) + u(i)*s(i-1))
t(i) = (r(i) - u(i)*t(i-1))/(v(i) + u(i)*s(i-1))
end do
y = t
do i = n-1, 1, -1
y(i) = s(i)*y(i+1) + t(i)
end do
end subroutine Tridiagonal_Solv
In that subroutine I want x and y to be outputs. I don't now how to do this.
I write a subroutine which has to solve differential equation numerically. Here it is
subroutine Tridiagonal_Solv(n, a, b, alpha,beta)
implicit none
integer :: n, i
double precision, intent(in) :: a, b, alpha, beta
real :: h
double precision, dimension(0:5) :: x, r, u, v, w, z, s, t, y
h = Step(n, a, b)
x(1) = a
r(1) = alpha
do i = 2, n-1
x(i) = a+i*h
r(i) = 2*sin(x(i))
end do
r = beta
x = b
do i = 2, n-1
u(i) = -1/(h*h)
end do
u = 2/h
v(1) = -3/(2*h)
do i = 2, n-1
v(i) = 2/(h*h) + 1
end do
v = 3/(2*h)
w(1) = 2/h
do i = 2, n-1
w(i) = -1/(h*h)
end do
s(0) = 0
t(0) = 0
do i = 1, n
s(i) = -w(i)/(v(i) + u(i)*s(i-1))
t(i) = (r(i) - u(i)*t(i-1))/(v(i) + u(i)*s(i-1))
end do
y = t
do i = n-1, 1, -1
y(i) = s(i)*y(i+1) + t(i)
end do
end subroutine Tridiagonal_Solv
In that subroutine I want x and y to be outputs. I don't now how to do this.