! Exact solution:
subroutine exact(t)
use work

implicit real(8) (a-h, o-z)
real(8) y(ngl), dy(ngl)
real(8) m_fi(ngl), fi_y(ngl), fi_dy(ngl)		! _ = *; fi(ngl, ngl)
real(8) mm, d, y_st
real(8) pp0(ngl)

real(8) S_1(ngl,ngl)

! Eigenvalues
! 1  0.9454183848254148    
! 2   3.478563691072702    
! 3  7.6017924101882628E-02

! Eigenvectors
!          1                       2                       3
!   1.000000000000000       1.000000000000000       1.000000000000000    
!  0.6273622311997391E-01  -1.237503200155199       1.920122969535155    
! -0.8665933062416533       1.064355928530766       1.292949362710746    

! Defineste Constante:
Const_If: IF(t ==tEnd) THEN

om(1) = 7.6017924101882628D-02
om(2) = 0.9454183848254148_8    
om(3) = 3.478563691072702_8    

om =dsqrt(om)


fi(1,:) =(/1,                     1,                      1/)
fi(2,:) =(/1.920122969535155_8,  0.6273622311997391D-01, -1.237503200155199_8/)
fi(3,:) =(/1.292949362710746_8, -0.8665933062416533_8,    1.064355928530766_8/)

! S_1:
!  0.707106781186547      -0.316227766016838       0.158113883008419     
!						   0.632455532033676      -0.316227766016838     
!												   0.790569415042095     

s_1(1,:) =(/  0.707106781186547_8,   -0.316227766016838_8,     0.158113883008419_8 /)
s_1(2,:) =(/  0._8,	  0.632455532033676_8,    -0.316227766016838_8 /)    
s_1(3,:) =(/  0._8,	  0._8,     0.790569415042095_8 /)


! Eigenvectors: FI =S_1*X
do i=1,ngl
fi(:,i) =matmul(S_1, fi(:,i))
enddo

!pp(1) =25_8
do igl =1, ngl
pp0(igl)=dot_product(fi(:, igl), pp)
enddo

aa =0.; m_fi =0.;

do igl =1,ngl
d =1./( 1-(omp/om(igl))**2 )

m_fi =matmul(M, fi(:, igl))
mm =dot_product(fi(:,igl), m_fi)

y_st =pp0(igl)/( mm*om(igl)**2)

aa(igl) =y_st*d
enddo


!*
write(10,*) 'Omega:', om

write(10, '(/,a)') ' FI:'
do i =1,ngl
write(10,*) fi(i,:)
enddo

write(10,'(/,a)') 'PP0:'
write(10,*) pp0

write(10,'(/,a)') ' A:'
write(10,*) aa
!*

ENDIF Const_If


! Solutia:
do igl =1,ngl
y(igl)  =aa(igl)*dsin(omp*t) -omp*aa(igl)/om(igl)*dsin(om(igl)*t)
dy(igl) =omp*aa(igl)*( dcos(omp*t) -dcos(om(igl)*t) )
enddo

uExact = matmul(fi, y)

uPrimeExact = matmul(fi, dy)


end subroutine exact
