! 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)



! Eigenvalues
! 1   917.0400956800341    
! 2   207.5974745328022    
! 3   1815.362429787164    

! Eigenvectors (of R)
!          1                       2                       3
!   1.000000000000000       1.000000000000000       1.000000000000000    
! -0.7472706718779475      0.9249016223752760      -2.864636719006213    
! -0.7668418199290550      0.4027533908560142       4.095576068071734    

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

om(1) = 207.5974745328029_8
om(2) = 917.0400956800345_8    
om(3) = 1815.362429787164_8    

om =dsqrt(om)


fi(1,:) =(/1,                     1,                     1/)
fi(2,:) =(/0.9249016224980822_8, -0.7472706738423111_8, -2.864636719006211_8/)
fi(3,:) =(/0.4027533907363419_8, -0.7668418154179999_8,  4.095576068071733_8/)

! Eigenvectors:
do i=1,ngl
fi(i,:) =fi(i,:)/dsqrt(M(i,i))
enddo


do i =1,ngl
fi(2:ngl,i) =fi(2:ngl,i)/fi(1,i)
enddo
fi(1,:) =1_8

!pp(1) =25_8
pp0 =matmul(transpose(fi), pp)

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
